      PROGRAM SPVDGX                                                                                                                
C     PROGRAM SPVDGX (SPU,        SPV,        SPX,          SPOUT,              E2                                                  
C    1                                        INPUT,        OUTPUT,     )       E2                                                  
C    2         TAPE11=SPU, TAPE12=SPV, TAPE13=SPX,   TAPE14=SPOUT,                                                                  
C    3                                 TAPE5 =INPUT, TAPE6 =OUTPUT)                                                                 
C     -------------------------------------------------------------             E2                                                  
C                                                                               E2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)       E2                                                  
C     NOV 06/95 - F.MAJAESS (REVISE DEFAULT PACKING DENSITY VALUE)                                                                  
C     FEB 15/94 - F.MAJAESS (REVISE FOR TRUNCATIONS > 99)                                                                           
C     JUL 22/92 - E. CHAN  (REPLACE UNFORMATTED I/O WITH I/O ROUTINES)                                                              
C     JUL 13/92 - E. CHAN  (DIMENSION SELECTED VARIABLES AS REAL*8)                                                                 
C     JAN 29/92 - E. CHAN  (CONVERT HOLLERITH LITERALS TO ASCII)                                                                    
C     AUG 13/85 - B.DUGAS. (AUGMENTER DIMENSION DE ALP, DALP, ETC.)                                                                 
C     NOV   /80 - S. LAMBERT                                                                                                        
C                                                                               E2                                                  
CSPVDGX  - SPECTRAL V DOT GRAD X                                        3  1 C GE1                                                  
C                                                                               E3                                                  
CAUTHOR  - S.LAMBERT                                                            E3                                                  
C                                                                               E3                                                  
CPURPOSE - COMPUTES THE SPHERICAL HARMONIC COEFFICIENTS OF A TERM OF            E3                                                  
C          THE FORM V DOT GRAD X FROM THE GLOBAL SPHERICAL HARMONIC             E3                                                  
C          COEFFICIENTS OF U, V, AND X.                                         E3                                                  
C          NOTE - MANY APPLICATIONS OF THIS PROGRAM WILL REQUIRE THE            E3                                                  
C                 DIVERGENCE TO BE REMOVED FROM THE WIND FIELD.                 E3                                                  
C                 MAXIMUM DIMENSIONS ARE $M$ WAVES, $IM1$ GAUSSIAN              E3                                                  
C                 LONGITUDES, AND $J$ GAUSSIAN LATITUDES.                       E3                                                  
C                                                                               E3                                                  
CINPUT FILES...                                                                 E3                                                  
C                                                                               E3                                                  
C      SPU   = GLOBAL SPECTRAL MODEL WIND U-COMPONENT                           E3                                                  
C      SPV   = GLOBAL SPECTRAL MODEL WIND V-COMPONENT                           E3                                                  
C      SPX   = GLOBAL SPECTRAL FIELDS OF SOME VARIABLE.                         E3                                                  
C                                                                               E3                                                  
COUTPUT FILE...                                                                 E3                                                  
C                                                                               E3                                                  
C      SPOUT = GLOBAL SPECTRAL COEFF OF V DOT GRAD X.                           E3                                                  
C                                                                                                                                   
CINPUT PARAMETERS...                                                                                                                
C                                                                               E5                                                  
C      NLAT  = NUMBER OF LATITUDES IN THE TRANSFORM GAUSSIAN GRID.              E5                                                  
C      ILONG = LENGTH OF A GAUSSIAN GRID ROW.                                   E5                                                  
C                                                                               E5                                                  
CEXAMPLE OF INPUT CARD...                                                       E5                                                  
C                                                                               E5                                                  
C*  SPVDGX   52   64                                                            E5                                                  
C----------------------------------------------------------------------------                                                       
C                                                                                                                                   
C     * ALL SPECTRAL INPUT ARRAYS ARE IN SP, THEIR ORDER BEING                                                                      
C     * U, V, X AND ANS IN A CONTIGUOUS MANNER.                                                                                     
C     * ALL SLICES ARRAYS ARE IN SLICE, THEIR ORDER BEING U, V,                                                                     
C     * XLAT, XLNG, X AND ANS IN A CONTIGUOUS MANNER.                                                                               
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      COMPLEX SP(583440),F(2210)                                                                                                    
      DIMENSION SLICE(194,396)                                                                                                      
      REAL*8    ALP(2276),DALP(2276),EPSI(2276)                                                                                     
      REAL*8    S(96),WEIGHT(96),SIA(96),RAD(96),WOCS(96)                                                                           
      DIMENSION TRIGS(194)                                                                                                          
      DIMENSION LSR(2,66), LRW(2,66),IB(8,66)                                                                                       
      DIMENSION WRKS(12416),IFAX(10)                                                                                                
C                                                                                                                                   
      COMMON/BUFCOM/IBUF(8),IDAT(4420)                                                                                              
      LOGICAL OK                                                                                                                    
      DATA MAXX,MAXLG,MAXSC,MAXTMS/4420,194,583440,66/                                                                              
C-----------------------------------------------------------------------                                                            
      NFF=6                                                                                                                         
      CALL JCLPNT(NFF,11,12,13,14,5,6)                                                                                              
      DO 70 N=11,14                                                                                                                 
         REWIND N                                                                                                                   
   70 CONTINUE                                                                                                                      
C                                                                                                                                   
C     * GET SIZE OF SPECTRAL FIELDS.                                                                                                
C                                                                                                                                   
      CALL FBUFFIN(13,IBUF,-8,K,LEN)                                                                                                
      IF (K.GE.0) GOTO 900                                                                                                          
      BACKSPACE 13                                                                                                                  
      BACKSPACE 13                                                                                                                  
      IF (IBUF (1).NE.NC4TO8("SPEC")) CALL         XIT('SPVDGX',-1)                                                                 
      LRLMT=IBUF(7)                                                                                                                 
      IF(LRLMT.LT.100) CALL FXLRLMT (LRLMT,IBUF(5),IBUF(6),0)                                                                       
      CALL DIMGT(LSR,LA,LR,LM,KTR,LRLMT)                                                                                            
      IR=LM-1                                                                                                                       
      LAW=LA+LM                                                                                                                     
C                                                                                                                                   
C     * LRW IS RELATIVE POSITION OF SPECTRAL WAVE NUMBER FOR                                                                        
C     * MODEL WINDS.                                                                                                                
C                                                                                                                                   
      ALP(LAW+1)=0.E0                                                                                                               
      LMP1=LM+1                                                                                                                     
      DO 130 M=1,LMP1                                                                                                               
         LRW(1,M)=LSR(2,M)                                                                                                          
         LRW(2,M)=LSR(2,M)                                                                                                          
  130 CONTINUE                                                                                                                      
C                                                                                                                                   
C     * READ SIZE OF TRANSFORM GRID                                                                                                 
C                                                                                                                                   
      READ(5,5000,END=910) ILAT,LON                                             E4                                                  
      IF (LON.GT.MAXLG-2) THEN                                                                                                      
         WRITE (6,6015) LON                                                                                                         
         LON=MAXLG-2                                                                                                                
      ENDIF                                                                                                                         
      WRITE(6,6010) LON,ILAT                                                                                                        
C                                                                                                                                   
C     * CALCULATE CONSTANTS.                                                                                                        
C                                                                                                                                   
      ILATH=ILAT/2                                                                                                                  
      ILH=MAXLG/2                                                                                                                   
      CALL GAUSSG(ILATH,S,WEIGHT,SIA,RAD,WOCS)                                                                                      
      CALL TRIGL (ILATH,S,WEIGHT,SIA,RAD,WOCS)                                                                                      
      CALL EPSCAL(EPSI,LSR,LM)                                                                                                      
      CALL FTSETUP(TRIGS,IFAX,LON)                                                                                                  
C                                                                                                                                   
      MAXLEV=MAXSC/(2*LAW+2*LA)                                                                                                     
      IF (MAXLEV.GT.MAXTMS) MAXLEV=MAXTMS                                                                                           
C---------------------------------------------------------------------                                                              
C                                                                                                                                   
C     * READ AS MANY FIELDS AS POSSIBLE.                                                                                            
C                                                                                                                                   
      NR=0                                                                                                                          
  150 ILEV=0                                                                                                                        
      NX=1                                                                                                                          
      DO 160 L=1,MAXLEV                                                                                                             
         CALL GETFLD2 (11,SP(NX),NC4TO8("SPEC"),-1,NC4TO8("   U"),                                                                  
     +                                            -1,IBUF,MAXX,OK)                                                                  
         IF(.NOT.OK)THEN                                                                                                            
           IF(ILEV+NR.EQ.0) CALL                   XIT('SPVDGX',-2)                                                                 
           IF(ILEV.EQ.0) THEN                                                                                                       
             WRITE (6,6020) NR                                                                                                      
             CALL                                  XIT('SPVDGX',0)                                                                  
           ENDIF                                                                                                                    
           GO TO 170                                                                                                                
         ENDIF                                                                                                                      
         IF(NR+L.EQ.1) THEN                                                                                                         
           CALL PRTLAB (IBUF)                                                                                                       
           NPACK=MIN(2,IBUF(8))                                                                                                     
         ENDIF                                                                                                                      
         NX=NX+LAW                                                                                                                  
         ILEV=ILEV+1                                                                                                                
         NR=NR+1                                                                                                                    
  160 CONTINUE                                                                                                                      
C                                                                                                                                   
  170 DO 180 L=1,ILEV                                                                                                               
         CALL GETFLD2 (12,SP(NX),NC4TO8("SPEC"),-1,NC4TO8("   V"),                                                                  
     +                                            -1,IBUF,MAXX,OK)                                                                  
         IF (.NOT.OK) CALL                         XIT('SPVDGX',-3)                                                                 
         IF (NR+L.EQ.ILEV+1) CALL PRTLAB (IBUF)                                                                                     
         IF(NR+L.EQ.1) NPACK=MIN(NPACK,IBUF(8))                                                                                     
         NX=NX+LAW                                                                                                                  
  180 CONTINUE                                                                                                                      
C                                                                                                                                   
      DO 190 L=1,ILEV                                                                                                               
         CALL GETFLD2 (13,SP(NX),NC4TO8("SPEC"),-1,0,-1,IBUF,MAXX,OK)                                                               
         IF (.NOT.OK) CALL                         XIT('SPVDGX',-4)                                                                 
         IF (NR+L.EQ.ILEV+1) CALL PRTLAB (IBUF)                                                                                     
         IF(NR+L.EQ.1) NPACK=MIN(NPACK,IBUF(8))                                                                                     
         NX=NX+LA                                                                                                                   
         DO 190 N=1,8                                                                                                               
            IB(N,L)=IBUF(N)                                                                                                         
  190 CONTINUE                                                                                                                      
C                                                                                                                                   
C     * DETERMINE POSITION OF SPECTRAL FIELDS.                                                                                      
C                                                                                                                                   
      ISU=1                                                                                                                         
      ISV=ISU+LAW*ILEV                                                                                                              
      ISX=ISV+LAW*ILEV                                                                                                              
      ISA=ISX+LA*ILEV                                                                                                               
C                                                                                                                                   
C     * DETERMINE POSITION OF FIELDS ON SLICE ARRAY                                                                                 
C                                                                                                                                   
      IGU=1                                                                                                                         
      IGV=IGU+ILEV                                                                                                                  
      ILT=IGV+ILEV                                                                                                                  
      ILG=ILT+ILEV                                                                                                                  
      IGX=ILG+ILEV                                                                                                                  
      IGA=IGX+ILEV                                                                                                                  
C                                                                                                                                   
C     * INITIALISE ANSWER TO ZERO AT START OF COMPUTATION.                                                                          
C                                                                                                                                   
      DO 200 J=ISA,ISA+LA*ILEV                                                                                                      
         SP(J)=CMPLX(0.E0,0.E0)                                                                                                     
  200 CONTINUE                                                                                                                      
C                                                                                                                                   
C      LOOP OVER TRANSFORM GRID LATITUDES.                                                                                          
C                                                                                                                                   
      DO 600 IH=1,ILAT                                                                                                              
         WLJ=WEIGHT(IH)                                                                                                             
         CALL ALPST2(ALP,LSR,LM,S(IH),EPSI)                                                                                         
         CALL ALPDY2(DALP,ALP,LSR,LM,EPSI)                                                                                          
C                                                                                                                                   
C     * COMPUTE FOURIER COEFFICIENTS AT ALL LEVELS FIRST FOR                                                                        
C     * U,V AND THEN FOR X AND XLNG.                                                                                                
C                                                                                                                                   
         NUMBER=2*ILEV                                                                                                              
         CALL STAF (SLICE(1,IGU),SP(ISU),LRW,LM,LAW,ILH,NUMBER,ALP)                                                                 
         CALL STAF (SLICE(1,IGX),SP(ISX),LSR,LM,LA,ILH,ILEV,ALP)                                                                    
         CALL STAF (SLICE(1,ILT),SP(ISX),LSR,LM,LA,ILH,ILEV,DALP)                                                                   
C                                                                                                                                   
         DO 520 M=1,LM                                                                                                              
            AM=FLOAT(M-1)                                                                                                           
            IREAL=2*M-1                                                                                                             
            IMAG=IREAL+1                                                                                                            
C                                                                                                                                   
C                                                                                                                                   
C                                                                                                                                   
            DO 520 L=0,ILEV-1                                                                                                       
              SLICE(IREAL,ILG+L)=-AM*SLICE(IMAG,IGX+L)                                                                              
              SLICE(IMAG,ILG+L) = AM*SLICE(IREAL,IGX+L)                                                                             
  520    CONTINUE                                                                                                                   
C                                                                                                                                   
C     * COMPUTE GRID POINT VALUES                                                                                                   
C                                                                                                                                   
         NUMBER=4*ILEV                                                                                                              
         CALL FFGFW(SLICE,MAXLG,SLICE,ILH,IR,LON,WRKS,NUMBER,                                                                       
     1              IFAX,TRIGS)                                                                                                     
         COSSQ=1.E0-S(IH)*S(IH)                                                                                                     
         DO 540 L=0,ILEV-1                                                                                                          
C                                                                                                                                   
C     * COMPUTE PRODUCT ON THE GRID                                                                                                 
C                                                                                                                                   
         DO 540 I=1,LON                                                                                                             
            A=SLICE(I,IGU+L)*SLICE(I,ILG+L)                                                                                         
            B=SLICE(I,IGV+L)*SLICE(I,ILT+L)                                                                                         
            SLICE(I,IGA+L)=(A+B)/COSSQ                                                                                              
  540    CONTINUE                                                                                                                   
C                                                                                                                                   
C     * COMPUTE FOURIER COEFFICIENTS                                                                                                
C                                                                                                                                   
         CALL FFWFG (SLICE(1,IGA),ILH,SLICE(1,IGA),MAXLG,IR,LON,WRKS,                                                               
     1               ILEV,IFAX,TRIGS)                                                                                               
C                                                                                                                                   
C     * COMPUTE SPHERICAL HARMONIC COEFFICIENTS                                                                                     
C                                                                                                                                   
         DO 550 L=1,LAW                                                                                                             
            ALP(L)=ALP(L)*WLJ                                                                                                       
  550    CONTINUE                                                                                                                   
         CALL FAST (SP(ISA),SLICE(1,IGA),LSR,LM,LA,ILH,ILEV,ALP)                                                                    
  600 CONTINUE                                                                                                                      
C                                                                                                                                   
C     * WRITE OUTPUT ON UNIT 14                                                                                                     
C                                                                                                                                   
      NX=ISA                                                                                                                        
      DO 700 L=1,ILEV                                                                                                               
       CALL SETLAB(IBUF,NC4TO8("SPEC"),IB(2,L),NC4TO8("VDGX"),                                                                      
     +                               IB(4,L),LA,1,LRLMT,NPACK)                                                                      
       CALL PUTFLD2 (14,SP(NX),IBUF,MAXX)                                                                                           
       NX=NX+LA                                                                                                                     
  700 CONTINUE                                                                                                                      
      GO TO 150                                                                                                                     
C                                                                                                                                   
C     * E.O.F. ON INPUT.                                                                                                            
C                                                                                                                                   
  900 CALL                                         XIT('SPVDGX',-5)                                                                 
  910 CALL                                         XIT('SPVDGX',-6)                                                                 
C-----------------------------------------------------------------------                                                            
 5000 FORMAT(10X,2I5)                                                           E4                                                  
 6010 FORMAT('0 LON,LAT=',2I6)                                                                                                      
 6015 FORMAT(' ** WARNING ** LON HAS BEEN REDUCED FROM',I5,                                                                         
     1       ' ** WARNING **')                                                                                                      
 6020 FORMAT('0SPVDGX PRODUCED',I6,' RECORDS')                                                                                      
      END                                                                                                                           
