      PROGRAM GWTQD                                                                                                                 
C     PROGRAM GWTQD (GGBIGU,       GGBIGV,       SPVORT,       SPDIV,           D2                                                  
C    1                                           INPUT,        OUTPUT,  )       D2                                                  
C    2         TAPE1=GGBIGU, TAPE2=GGBIGV, TAPE3=SPVORT, TAPE4=SPDIV,                                                               
C    3                                     TAPE5=INPUT,  TAPE6=OUTPUT)                                                              
C     ----------------------------------------------------------------          D2                                                  
C                                                                               D2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)       D2                                                  
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     FEB 22/85 - B.DUGAS. (SET TO 22 T20 GRIDS AT A TIME INSTEAD OF 25)                                                            
C                                                                               D2                                                  
CGWTQD   - COMPUTES SPECTRAL Q,D FROM GRID WINDS                        2  2 C  D1                                                  
C                                                                               D3                                                  
CAUTHOR  - J.D.HENDERSON                                                        D3                                                  
C                                                                               D3                                                  
CPURPOSE - COMPUTES THE SPECTRAL COEEFICIENT FILES FOR VORTICITY AND            D3                                                  
C          DIVERGENCE FROM WIND COMPONENTS ON GAUSSIAN GRIDS.                   D3                                                  
C          NOTE - WIND INPUT MAY BE REAL WIND OR MODEL WIND WHERE THE           D3                                                  
C                 MODEL WIND COMPONENTS = (U,V)*COS(LAT)/(EARTH RADIUS).        D3                                                  
C                                                                               D3                                                  
CINPUT FILES...                                                                 D3                                                  
C                                                                               D3                                                  
C      GGBIGU = GLOBAL GAUSSIAN GRIDS OF U WIND COMPONENTS.                     D3                                                  
C      GGBIGV = GLOBAL GAUSSIAN GRIDS OF V WIND COMPONENTS.                     D3                                                  
C                                                                               D3                                                  
COUTPUT FILES...                                                                D3                                                  
C                                                                               D3                                                  
C      SPVORT = SPECTRAL VORTICITY                                              D3                                                  
C      SPDIV  = SPECTRAL DIVERGENCE                                             D3                                                  
C                                                                                                                                   
CINPUT PARAMETERS...                                                                                                                
C                                                                               D5                                                  
C      LRT,LMT = N,M TRUNCATION WAVE NUMBERS                                    D5                                                  
C      KTR     = SPECTRAL TRUNCATION TYPE                                       D5                                                  
C      KUV     = 0, ASSUMES MODEL WINDS,                                        D5                                                  
C              = 1, ASSUMES TRUE WINDS.                                         D5                                                  
C                                                                               D5                                                  
CEXAMPLE OF INPUT CARD...                                                       D5                                                  
C                                                                               D5                                                  
C*   GWTQD   20   20    2    0                                                  D5                                                  
C--------------------------------------------------------------------------                                                         
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      COMPLEX DQ(429000)                                                                                                            
      DIMENSION U(1852800),V(1852800),WRKL(39200),WRKS(39200)                                                                       
      DIMENSION IFAX(10),TRIGS(194)                                                                                                 
      REAL*8 SL(96),CL(96),WL(96),WOSSL(96),RAD(96)                                                                                 
      REAL*8 ALP(2275),DALP(2275),EPSI(2275)                                                                                        
C                                                                                                                                   
      LOGICAL OK                                                                                                                    
      INTEGER LSR(2,66),IBU(8,100),IBV(8,100)                                                                                       
      COMMON/ICOM/IBUF(8),IDAT(18528)                                                                                               
      DATA MAXX,MAXLG,MAXGG,MAXSP/18528,194,1852800,429000/                                                                         
      DATA MAXL/100/                                                                                                                
C-----------------------------------------------------------------------                                                            
      NFF=6                                                                                                                         
      CALL JCLPNT(NFF,1,2,3,4,5,6)                                                                                                  
      DO 50 N=1,4                                                                                                                   
   50 REWIND N                                                                                                                      
C                                                                                                                                   
C     * READ SPECTRAL ARRAY SIZE, TRUNCATION TYPE, AND KUV.                                                                         
C     * KUV=1 CONVERTS REAL WINDS TO MODEL WINDS.                                                                                   
C                                                                                                                                   
      READ(5,5010,END=904) LRT,LMT,KTR,KUV                                      D4                                                  
      IF(KUV.EQ.1) WRITE(6,6008)                                                                                                    
C     LRLMT=1000*(LRT+1)+10*(LMT+1)+KTR                                                                                             
      CALL FXLRLMT (LRLMT,LRT+1,LMT+1,KTR)                                                                                          
      CALL DIMGT(LSR,LA,LR,LM,KTR,LRLMT)                                                                                            
      CALL PRLRLMT (LA,LR,LM,KTR,LRLMT)                                                                                             
C                                                                                                                                   
C     * READ SIZE OF GRIDS AND DETERMINE MAXLEV.                                                                                    
C                                                                                                                                   
      CALL FBUFFIN(1,IBUF,-8,K,LEN)                                                                                                 
      BACKSPACE 1                                                                                                                   
      BACKSPACE 1                                                                                                                   
      IF (IBUF(1).NE.NC4TO8("GRID")) CALL          XIT('GWTQD',-1)                                                                  
      ILG1=IBUF(5)                                                                                                                  
      ILG=ILG1-1                                                                                                                    
      ILH=(ILG1+1)/2                                                                                                                
      ILAT=IBUF(6)                                                                                                                  
      NWDS=ILG1*ILAT                                                                                                                
      IF (ILG.GT.MAXLG-2) CALL                     XIT('GWTQD',-2)                                                                  
C                                                                                                                                   
      MAX1=MAXGG/NWDS                                                                                                               
      MAX2=MAXSP/(2*LA)                                                                                                             
      MAXLEV=MAX1                                                                                                                   
      IF (MAXLEV.GT.MAX2) MAXLEV=MAX2                                                                                               
      IF (MAXLEV.GT.MAXL) MAXLEV=MAXL                                                                                               
C                                                                                                                                   
C     * CALCULATE THE CONSTANTS.  GAUSSG COMPUTES GAUSSIAN LATITUDES,                                                               
C     * ETC.  TRIGL MAKES THEM GLOBAL. FTSETUP CALCULATES CONSTANTS FOR                                                             
C     * ECMWF VFFT'S.                                                                                                               
C                                                                                                                                   
      ILATH=ILAT/2                                                                                                                  
      CALL EPSCAL(EPSI,LSR,LM)                                                                                                      
      CALL GAUSSG(ILATH,SL,WL,CL,RAD,WOSSL)                                                                                         
      CALL TRIGL(ILATH,SL,WL,CL,RAD,WOSSL)                                                                                          
      CALL FTSETUP (TRIGS,IFAX,ILG)                                                                                                 
C-----------------------------------------------------------------------                                                            
C                                                                                                                                   
C     * READ AS MANY U AND V FIELDS AS POSSIBLE.                                                                                    
C                                                                                                                                   
      NR=0                                                                                                                          
      NPACK=2                                                                                                                       
  140 ILEV=0                                                                                                                        
      NX=1                                                                                                                          
      DO 160 L=1,MAXLEV                                                                                                             
         CALL GETFLD2 (1,U(NX),NC4TO8("GRID"),-1,-1,-1,IBUF,MAXX,OK)                                                                
         IF (.NOT.OK)THEN                                                                                                           
           IF (ILEV+NR.EQ.0)THEN                                                                                                    
             CALL                                  XIT('GWTQD',-3)                                                                  
           ELSE                                                                                                                     
             IF (ILEV.EQ.0)THEN                                                                                                     
               WRITE (6,6010) NR                                                                                                    
               CALL                                XIT('GWTQD',0)                                                                   
             ELSE                                                                                                                   
               GO TO 170                                                                                                            
             ENDIF                                                                                                                  
           ENDIF                                                                                                                    
         ENDIF                                                                                                                      
         IF(NR.EQ.0) THEN                                                                                                           
           NPACK=MIN(NPACK,IBUF(8))                                                                                                 
         ENDIF                                                                                                                      
         IF (NR+L.EQ.1)                CALL PRTLAB (IBUF)                                                                           
                                       NX=NX+NWDS                                                                                   
                                       ILEV=ILEV+1                                                                                  
                                       NR=NR+1                                                                                      
         DO 160 N=1,8                                                                                                               
            IBU(N,L)=IBUF(N)                                                                                                        
  160 CONTINUE                                                                                                                      
C                                                                                                                                   
  170 NX=1                                                                                                                          
      DO 180 L=1,ILEV                                                                                                               
         CALL GETFLD2 (2,V(NX),NC4TO8("GRID"),-1,-1,-1,IBUF,MAXX,OK)                                                                
         IF (.NOT.OK) CALL                         XIT('GWTQD',-4)                                                                  
         IF (NR+L.EQ.ILEV+1)           CALL PRTLAB (IBUF)                                                                           
                                       NX=NX+NWDS                                                                                   
         IF(NR.EQ.0) THEN                                                                                                           
           NPACK=MIN(NPACK,IBUF(8))                                                                                                 
         ENDIF                                                                                                                      
         DO 180 N=1,8                                                                                                               
            IBV(N,L)=IBUF(N)                                                                                                        
  180 CONTINUE                                                                                                                      
C                                                                                                                                   
      DO 190 L=1,ILEV                                                                                                               
         CALL CMPLBL(0,IBU(1,L),0,IBV(1,L),OK)                                                                                      
         IF(.NOT.OK)THEN                                                                                                            
            CALL PRTLAB (IBU(1,L))                                                                                                  
            CALL PRTLAB (IBV(1,L))                                                                                                  
           CALL                                    XIT('GWTQD',-10-L)                                                               
         ENDIF                                                                                                                      
  190 CONTINUE                                                                                                                      
C                                                                                                                                   
C     * COMPUTE VORT AND DIV FROM WINDS.                                                                                            
C                                                                                                                                   
      IF(KUV.EQ.1)THEN                                                                                                              
        CALL LWBW (U,ILG1,ILAT,ILEV,CL,1)                                                                                           
        CALL LWBW (V,ILG1,ILAT,ILEV,CL,1)                                                                                           
      ENDIF                                                                                                                         
      CALL GWAQD (DQ,LSR,LA,LM,U,V,ILG1,ILH,ILAT,SL,WOSSL,ALP,                                                                      
     1            DALP,EPSI,WRKL,WRKS,ILEV,IFAX,TRIGS)                                                                              
C                                                                                                                                   
C     * PUT SPECTRAL VORT,DIV ON FILES 3,4.                                                                                         
C                                                                                                                                   
      ND=1                                                                                                                          
      NQ=LA*ILEV+1                                                                                                                  
      DO 200 L=1,ILEV                                                                                                               
         ITIM=IBU(2,L)                                                                                                              
         LEVL=IBU(4,L)                                                                                                              
         CALL SETLAB(IBUF,NC4TO8("SPEC"),ITIM,NC4TO8("VORT"),LEVL,                                                                  
     +                                           LA,1,LRLMT,NPACK)                                                                  
         CALL PUTFLD2(3,DQ(NQ),IBUF,MAXX)                                                                                           
         IF(NR+L.EQ.ILEV+1) CALL PRTLAB (IBUF)                                                                                      
                                                                                                                                    
C                                                                                                                                   
         IBUF(3)=NC4TO8(" DIV")                                                                                                     
         CALL PUTFLD2(4,DQ(ND),IBUF,MAXX)                                                                                           
         IF(NR+L.EQ.ILEV+1) CALL PRTLAB (IBUF)                                                                                      
C                                                                                                                                   
         NQ=NQ+LA                                                                                                                   
         ND=ND+LA                                                                                                                   
C                                                                                                                                   
  200 CONTINUE                                                                                                                      
      GO TO 140                                                                                                                     
C                                                                                                                                   
C     * E.O.F. ON INPUT.                                                                                                            
C                                                                                                                                   
  904 CALL                                         XIT('GWTQD',-5)                                                                  
C-----------------------------------------------------------------------                                                            
 5010 FORMAT(10X,4I5)                                                           D4                                                  
 6008 FORMAT('0 INCLUDE WIND CONVERSION'/)                                                                                          
 6010 FORMAT('0GWTQD CONVERTED',I6,' RECORDS')                                                                                      
      END                                                                                                                           
