      PROGRAM GGACOF                                                                                                                
C     PROGRAM GGACOF (GGFILE,       SPFILE,       INPUT,       OUTPUT,  )       D2                                                  
C    1          TAPE1=GGFILE, TAPE2=SPFILE, TAPE5=INPUT, TAPE6=OUTPUT)                                                              
C     ----------------------------------------------------------------          D2                                                  
C                                                                               D2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)       D2                                                  
C     FEB 15/94 - F.MAJAESS (REVISE FOR TRUNCATIONS > 99)                                                                           
C     JUL 13/92 - E. CHAN  (DIMENSION SELECTED VARIABLES AS REAL*8)                                                                 
C     JAN 29/92 - E. CHAN  (CONVERT HOLLERITH LITERALS TO ASCII)                                                                    
C     NOV 21/90 - F.MAJAESS. (SKIP CARDS OTHER THAN "GRID" AND "ZONL")                                                              
C     DEC 17/87 - F.MAJAESS. (FIX THE PROGRAM SO THAT IT WON'T ABORT FOR                                                            
C                             THE ZONAL CASE WITH M DIFFERENT FROM ZERO).                                                           
C     AUG 13/85 - B.DUGAS.   (AUGMENTER DIM SPECTRALE A T42)                                                                        
C     MAR 31/81 - J.D.HENDERSON                                                                                                     
C                                                                               D2                                                  
CGGACOF  - CONVERTS GAUSSIAN GRID (OR ZONAL CROSS-SECTION) FILE TO              D1                                                  
C          SPECTRAL COEFFICIENT FILE                                    1  1 C GD1                                                  
C                                                                               D3                                                  
CAUTHOR  - J.D.HENDERSON                                                        D3                                                  
C                                                                               D3                                                  
CPURPOSE - CONVERTS A FILE OF GLOBAL GAUSSIAN GRIDS (OR GAUSSIAN ZONAL          D3                                                  
C          CROSS-SECTIONS) TO GLOBAL SPECTRAL COEFFICIENTS PACKED 2 TO 1.       D3                                                  
C          NOTE - THE ZONAL CASE OUTPUT SPECTRAL ARRAY HAS NON-ZERO VALUES      D3                                                  
C                 ONLY FOR THE M=0 PART (ESSENTIALY, THE LEGENDRE TRANSFORM).   D3                                                  
C                 ALL SPECTRAL FIELDS WILL BE OF THE SAME SIZE.                 D3                                                  
C                                                                               D3                                                  
CINPUT FILE...                                                                  D3                                                  
C                                                                               D3                                                  
C      GGFILE = FILE OF GAUSSIAN GRIDS (GLOBAL OR ZONAL)                        D3                                                  
C                                                                               D3                                                  
COUTPUT FILE...                                                                 D3                                                  
C                                                                               D3                                                  
C      SPFILE = GLOBAL SPECTRAL COEFFICIENTS.                                   D3                                                  
C                                                                                                                                   
CINPUT PARAMETERS...                                                                                                                
C                                                                               D5                                                  
C      LRT,LMT = N,M TRUNCATION SPECTRAL WAVE NUMBERS                           D5                                                  
C      KTR     = TRUNCATION TYPE (0=RHOMBOIDAL, 2=TRIANGULAR)                   D5                                                  
C      KUV     = 0 FOR NORMAL ANALYSIS                                          D5                                                  
C                1 TO CONVERT REAL WINDS TO MODEL WINDS                         D5                                                  
C      KPACK   = SPECTRAL PACKING DENSITY (0 DEFAULTS TO 2)                     D5                                                  
C                                                                               D5                                                  
CEXAMPLE OF INPUT CARD...                                                       D5                                                  
C                                                                               D5                                                  
C*  GGACOF   20   20    2    1    2                                             D5                                                  
C---------------------------------------------------------------------------                                                        
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      LOGICAL OK                                                                                                                    
      INTEGER IB(8,100)                                                                                                             
      INTEGER LSR(2,66),IFAX(10)                                                                                                    
      REAL TMPA(32016),WRKS(12416),WRKL(19600),TRIGS(194)                                                                           
      REAL*8 SL(96),CL(96),WL(96),WOSSL(96),RAD(96)                                                                                 
      REAL*8 ALP(2275),EPSI(2275)                                                                                                   
                                                                                                                                    
      COMMON/W/W(2304400)                                                                                                           
      COMMON/ICOM/IBUF(8),IDAT(18528)                                                                                               
                                                                                                                                    
      DATA MAXX,MAXW,MAXLG,MAXT/18528,2304400,194,32016/                                                                            
      DATA MAXL/100/                                                                                                                
      EQUIVALENCE (TMPA(1),WRKS(1)),(TMPA(12417),WRKL(1))                                                                           
C                                                                                                                                   
C---------------------------------------------------------------------                                                              
      NF=4                                                                                                                          
      CALL JCLPNT(NF,1,2,5,6)                                                                                                       
      REWIND 1                                                                                                                      
      REWIND 2                                                                                                                      
                                                                                                                                    
C     * READ SPECTRAL N AND M TRUNCATION NUMBERS.                                                                                   
C     * NORMALLY USE (IR,IR). FOR WINDS USE (IR+1,IR).                                                                              
C     * KUV=1 CONVERTS REAL WINDS TO MODEL WINDS.                                                                                   
C     * KPACK = OUTPUT PACKING DENSITY (0 DEFAULTS TO 2).                                                                           
                                                                                                                                    
      READ(5,5010,END=901) LRT,LMT,KTR,KUV,KPACK                                D4                                                  
      IF(KPACK.EQ.0) KPACK=2                                                                                                        
      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     * READ THE FIRST GRID FIELD TO GET THE SIZE.                                                                                  
                                                                                                                                    
      NR=0                                                                                                                          
  100 CALL GETFLD2(1,W( 1 ),-1,-1,-1,-1,IBUF,MAXX,OK)                                                                               
      IF(.NOT.OK)THEN                                                                                                               
        WRITE(6,6010) NR                                                                                                            
        CALL                                       XIT('GGACOF',-1)                                                                 
      ENDIF                                                                                                                         
      IF (IBUF(1).NE.NC4TO8("GRID") .AND.                                                                                           
     +    IBUF(1).NE.NC4TO8("ZONL")      ) GO TO 100                                                                                
      BACKSPACE  1                                                                                                                  
      BACKSPACE  1                                                                                                                  
      IF (IBUF(1).EQ.NC4TO8("GRID"))   THEN                                                                                         
                                                                                                                                    
C**********************************************************************                                                             
C         * GRID FIELD CASE.                                                                                                        
                                                                                                                                    
          ILG1=IBUF(5)                                                                                                              
          IF (ILG1+1.GT.MAXLG)   THEN                                                                                               
              WRITE(6,6005) ILG1,MAXLG-1                                                                                            
              CALL                                 XIT('GGACOF',-2)                                                                 
          ENDIF                                                                                                                     
          MAXLG=ILG1+1                                                                                                              
          ILG=ILG1-1                                                                                                                
          ILAT=IBUF(6)                                                                                                              
          LGG=ILG1*ILAT                                                                                                             
          ILATH=ILAT/2                                                                                                              
          CALL PRTLAB (IBUF)                                                                                                        
                                                                                                                                    
C         * FIRST TIME ONLY, SET CONSTANTS FOR SPECTRAL-GRID TRANSFOR-                                                              
C         * MATION. GAUSSG COMPUTES THE VALUE OF THE GAUSSIAN LATITUDES                                                             
C         * AND THEIR SINES AND COSINES. TRIGL MAKES THEM GLOBAL (S TO N).                                                          
                                                                                                                                    
          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         * FIND HOW MANY SP,GG PAIRS WE CAN FIT (MAXL).                                                                            
                                                                                                                                    
          NWDS=2*LA+LGG                                                                                                             
          MAXLEV=MAXW/NWDS                                                                                                          
          IF(MAXLEV.GT.MAXL) MAXLEV=MAXL                                                                                            
          NGG1=2*LA*MAXLEV+1                                                                                                        
          WRITE(6,6015) LR,LM,ILG1,ILAT,MAXLEV                                                                                      
                                                                                                                                    
C--------------------------------------------------------------------------                                                         
C         * READ AS MANY GRIDS AS POSSIBLE.                                                                                         
C         * KUV=1 CONVERTS REAL WINDS TO MODEL WINDS.                                                                               
                                                                                                                                    
  110     ILEV=0                                                                                                                    
          NGG=NGG1                                                                                                                  
          DO 160 L=1,MAXLEV                                                                                                         
          CALL GETFLD2(1,W(NGG),NC4TO8("GRID"),-1,-1,-1,IBUF,MAXX,OK)                                                               
          IF(.NOT.OK)THEN                                                                                                           
            IF(ILEV+NR.EQ.0) CALL                  XIT('GGACOF',-3)                                                                 
            IF(ILEV   .EQ.0) GO TO 510                                                                                              
            GO TO 180                                                                                                               
          ENDIF                                                                                                                     
          DO 150 I=1,8                                                                                                              
  150     IB(I,L)=IBUF(I)                                                                                                           
          ILEV=ILEV+1                                                                                                               
          NGG=NGG+LGG                                                                                                               
  160     CONTINUE                                                                                                                  
                                                                                                                                    
C         * CONVERT ILEV GRIDS TO SPECTRAL COEFF.                                                                                   
C         * NORMAL SPEC FIELDS ARE (LR,LM) BUT WINDS ARE (LRW,LM).                                                                  
                                                                                                                                    
  180     IF(KUV.EQ.1) CALL LWBW (W(NGG1),ILG1,ILAT,ILEV,CL,1)                                                                      
          CALL GGAST (W(1),LSR,LM,LA, W(NGG1),ILG1,ILAT,ILEV,SL,WL,                                                                 
     1                           ALP,EPSI,WRKS,WRKL,TRIGS,IFAX,MAXLG)                                                               
                                                                                                                                    
C         * WRITE ALL THE SPECTRAL FIELDS ONTO FILE 2.                                                                              
                                                                                                                                    
          NSP=1                                                                                                                     
          DO 490 L=1,ILEV                                                                                                           
             CALL SETLAB(IBUF,NC4TO8("SPEC"),IB(2,L),IB(3,L),                                                                       
     +                              IB(4,L),LA,1,LRLMT,KPACK)                                                                       
             CALL PUTFLD2(2,W(NSP),IBUF,MAXX)                                                                                       
             NSP=NSP+2*LA                                                                                                           
             NR=NR+1                                                                                                                
  490     CONTINUE                                                                                                                  
                                                                                                                                    
C         * STOP IF ALL GRIDS HAVE BEEN PROCESSED.                                                                                  
C         * OTHERWISE GO BACK FOR THE NEXT SET.                                                                                     
                                                                                                                                    
  510     IF(.NOT.OK)THEN                                                                                                           
            WRITE(6,6010) NR                                                                                                        
            CALL                                   XIT('GGACOF',0)                                                                  
          ENDIF                                                                                                                     
          GO TO 110                                                                                                                 
                                                                                                                                    
      ELSE IF (IBUF(1).EQ.NC4TO8("ZONL"))   THEN                                                                                    
                                                                                                                                    
C***********************************************************************                                                            
C         * ZONAL CROSS-SECTION CASE.                                                                                               
                                                                                                                                    
          LAZ=LA                                                                                                                    
          LMZ=LM                                                                                                                    
          LRLMTZ=LRLMT                                                                                                              
          IF(LMT.NE.0)THEN                                                                                                          
            LMT=0                                                                                                                   
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)                                                                                      
          ENDIF                                                                                                                     
                                                                                                                                    
          ILAT=IBUF(5)                                                                                                              
          ILATH=ILAT/2                                                                                                              
          CALL PRTLAB (IBUF)                                                                                                        
                                                                                                                                    
C         * FIRST TIME ONLY, SET CONSTANTS FOR SPECTRAL-GRID TRANSFOR-                                                              
C         * MATION. GAUSSG COMPUTES THE VALUE OF THE GAUSSIAN LATITUDES                                                             
C         * AND THEIR SINES AND COSINES. TRIGL MAKES THEM GLOBAL (S TO N).                                                          
                                                                                                                                    
          CALL EPSCAL(EPSI,LSR,LM)                                                                                                  
          CALL GAUSSG(ILATH,SL,WL,CL,RAD,WOSSL)                                                                                     
          CALL  TRIGL(ILATH,SL,WL,CL,RAD,WOSSL)                                                                                     
                                                                                                                                    
C         * FIND HOW MANY SP,ZX PAIRS WE CAN FIT (MAXL).                                                                            
                                                                                                                                    
          NWDS=2*LA+ILAT                                                                                                            
          MAXLEV=MAXW/NWDS                                                                                                          
          IF(MAXLEV.GT.MAXL) MAXLEV=MAXL                                                                                            
          NZX1=2*LA*MAXLEV+1                                                                                                        
          WRITE(6,6115) LR,LMZ,ILAT,MAXLEV                                                                                          
                                                                                                                                    
C--------------------------------------------------------------------------                                                         
C         * READ AS MANY CROSS-SECTIONS AS POSSIBLE.                                                                                
C         * KUV=1 CONVERTS REAL WINDS TO MODEL WINDS.                                                                               
                                                                                                                                    
  610     ILEV=0                                                                                                                    
          NZX=NZX1                                                                                                                  
          DO 660 L=1,MAXLEV                                                                                                         
          CALL GETFLD2(1,W(NZX),NC4TO8("ZONL"),-1,-1,-1,IBUF,MAXX,OK)                                                               
          IF(.NOT.OK)THEN                                                                                                           
            IF(ILEV+NR.EQ.0) CALL                  XIT('GGACOF',-4)                                                                 
            IF(ILEV   .EQ.0) GO TO 710                                                                                              
            GO TO 680                                                                                                               
          ENDIF                                                                                                                     
          DO 650 I=1,8                                                                                                              
  650     IB(I,L)=IBUF(I)                                                                                                           
          ILEV=ILEV+1                                                                                                               
          NZX=NZX+ILAT                                                                                                              
  660     CONTINUE                                                                                                                  
                                                                                                                                    
C         * CONVERT ILEV FILEDS TO SPECTRAL COEFF.                                                                                  
C         * NORMAL SPEC FIELDS ARE (LR,LM) BUT WINDS ARE (LRW,LM).                                                                  
                                                                                                                                    
  680     IF (KUV.EQ.1)    THEN                                                                                                     
              A     =6.371E6                                                                                                        
              ISTART=NZX1-ILAT-1                                                                                                    
              DO 690 L=1,ILEV                                                                                                       
                 ISTART=ISTART+ILAT                                                                                                 
                 DO 690 J=1,ILAT                                                                                                    
                    W(J+ISTART)=W(J+ISTART)*CL(J)/A                                                                                 
  690         CONTINUE                                                                                                              
          ENDIF                                                                                                                     
          CALL ZXAST (W(1),LSR,LA, W(NZX1),ILAT,ILEV,SL,WL,ALP,EPSI)                                                                
                                                                                                                                    
C         * ABORT IF THE SIZE OF THE TEMPORARY ARRAY 'TMPA' IS NOT                                                                  
C         * SUFFICIENT.                                                                                                             
                                                                                                                                    
          IF(2*LAZ.GT.MAXT) CALL                   XIT('GGACOF',-5)                                                                 
                                                                                                                                    
C         * SET THE NEEDED ELEMENTS OF THE TEMPORARY ARRAY 'TMPA' TO ZERO.                                                          
                                                                                                                                    
          DO 695 I=1,2*LAZ                                                                                                          
             TMPA(I)=0.0E0                                                                                                          
  695     CONTINUE                                                                                                                  
                                                                                                                                    
C         * WRITE ALL THE SPECTRAL FIELDS ONTO FILE 2.                                                                              
                                                                                                                                    
          NSP=1                                                                                                                     
          DO 705 L=1,ILEV                                                                                                           
                                                                                                                                    
             CALL SETLAB(IBUF,NC4TO8("SPEC"),IB(2,L),IB(3,L),                                                                       
     +                            IB(4,L),LAZ,1,LRLMTZ,KPACK)                                                                       
                                                                                                                                    
C         * TRANSFER THE DATA TO BE WRITTEN INTO THE TEMPORARY ARRAY.                                                               
C         * THEN WRITE THEM OUT.                                                                                                    
                                                                                                                                    
             DO 700 I=1,(2*LA)                                                                                                      
                TMPA(I)=W(NSP+I-1)                                                                                                  
  700        CONTINUE                                                                                                               
             CALL PUTFLD2(2,TMPA(1),IBUF,MAXX)                                                                                      
             NSP=NSP+2*LA                                                                                                           
             NR=NR+1                                                                                                                
  705     CONTINUE                                                                                                                  
                                                                                                                                    
C         * STOP IF ALL GRIDS HAVE BEEN PROCESSED.                                                                                  
C         * OTHERWISE GO BACK FOR THE NEXT SET.                                                                                     
                                                                                                                                    
  710     IF(.NOT.OK)THEN                                                                                                           
            WRITE(6,6010) NR                                                                                                        
            CALL                                   XIT('GGACOF',0)                                                                  
          ENDIF                                                                                                                     
          GO TO 610                                                                                                                 
                                                                                                                                    
      ENDIF                                                                                                                         
                                                                                                                                    
C***********************************************************************                                                            
C     * E.O.F. ON INPUT.                                                                                                            
                                                                                                                                    
  901 CALL                                         XIT('GGACOF',-6)                                                                 
C---------------------------------------------------------------------                                                              
 5010 FORMAT(10X,5I5)                                                           D4                                                  
 6005 FORMAT(' GRIDS TOO BIG,  ILG1, MAXLG = ',I5,', ',I5)                                                                          
 6008 FORMAT('0 INCLUDE WIND CONVERSION'/)                                                                                          
 6010 FORMAT('0',I6,'  RECORDS CONVERTED')                                                                                          
 6015 FORMAT(' SPEC=',2I5,'  GRID=',2I5,'  MAXLEV=',I5)                                                                             
 6115 FORMAT(' SPEC=',2I5,'  ZONL=',I5,'  MAXLEV=',I5)                                                                              
      END                                                                                                                           
