      PROGRAM GGALL                                                                                                                 
C     PROGRAM GGALL (GG,       VG,       LL,       VL,       FL,                D2                                                  
C                                               INPUT,       OUTPUT,    )       D2                                                  
C    1        TAPE11=GG,TAPE12=VG,TAPE13=LL,TAPE14=VL,TAPE15=FL,                                                                    
C                                         TAPE5=INPUT, TAPE6=OUTPUT)                                                                
C     --------------------------------------------------------------            D2                                                  
C                                                                               D2                                                  
C     MAR 25/13 - S.KHARIN (ADD AREA-WEIGHTED INTERPOLATION.                    D2                                                  
C                           IMPROVE INTERPOLATION NEAR THE POLES.               D2                                                  
C                           ALLOW ZONALLY AVERAGED INPUT FIELDS.                D2                                                  
C                           CORRECT TREATMENT OF VECTOR FIELDS.                 D2                                                  
C                           ADD OPTION FOR OUTPUT GAUSSIAN GRIDS.               D2                                                  
C                           ADD OPTION FOR PROPAGATING SUPERLABELS.             D2                                                  
C                           USE PACKING DENSITY OF INPUT FIELD AS DEFAULT.      D2                                                  
C                           USE BILINEAR INTERPOLATION AS DEFAULT.              D2                                                  
C                           TREAT SPECIAL VALUES 1E38 AS MISSING.               D2                                                  
C                           USE MORE ACCURATE PI. HARMONIZE CODE WITH LLAGG).   D2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)                                                           
C     FEB 09/94 - D. RAMSDEN (HANDLE CASES WITH SHIFTED LAT/LON GRIDS)                                                              
C     JAN 29/92 - E. CHAN  (CONVERT HOLLERITH LITERALS TO ASCII)                                                                    
C     MAY 11/83 - R.LAPRISE.                                                                                                        
C     MAR 23/81 - J.D.HENDERSON                                                                                                     
C                                                                               D2                                                  
CGGALL   - CONVERT GAUSSIAN GRID/ZONL TO LAT-LON OR GAUSSIAN GRID/ZONL  2  2 C GD1                                                  
C                                                                               D3                                                  
CAUTHOR  - J.D.HENDERSON                                                        D3                                                  
C                                                                               D3                                                  
CPURPOSE - CONVERTS GLOBAL GAUSSIAN GRIDS OR ZONAL AVERAGES TO GLOBAL LAT-LON   D3                                                  
C          OR GAUSSIAN GRIDS OR ZONAL AVERAGES BY BILINEAR, CUBIC, OR           D3                                                  
C          AREA-WEIGTED INTERPOLATION.                                          D3                                                  
C          NOTE - OUTPUT LAT-LON GRIDS MAY BE SHIFTED AND UNSHIFTED.            D3                                                  
C                 LAT-UNSHIFTED GRID INCLUDES THE POLES.                        D3                                                  
C                 LON-UNSHIFTED GRID INCLUDES AN EXTRA CYCLIC LONGITUDE.        D3                                                  
C                 SPVAL=1.E38 VALUES ARE TREATED AS MISSING.                    D3                                                  
C                                                                               D3                                                  
CINPUT FILES...                                                                 D3                                                  
C                                                                               D3                                                  
C      GG = GLOBAL GAUSSIAN GRIDS OR ZONAL AVERAGES OF A SCALAR FIELD           D3                                                  
C           OR U- OR V-COMPONENT OF A VECTOR FIELD.                             D3                                                  
C      VG = (OPTIONAL) GLOBAL GAUSSIAN GRIDS OR ZONAL AVERAGES OF THE           D3                                                  
C           V-COMPONENT OF A VECTOR FIELD (SEE INPUT CARD PARAMETER IFLD=2).    D3                                                  
C                                                                               D3                                                  
COUTPUT FILES...                                                                D3                                                  
C                                                                               D3                                                  
C      LL = GLOBAL LAT-LON OR GAUSSIAN GRIDS OR ZONAL AVERAGES OF A SCALAR      D3                                                  
C           FIELD OR U- OR V-COMPONENT OF A VECTOR FIELD.                       D3                                                  
C      VL = (OPTIONAL) GLOBAL LAT-LON OR GAUSSIAN GRIDS OR ZONAL AVERAGES OF    D3                                                  
C           THE V-COMPONENT OF A VECTOR FIELD.                                  D3                                                  
C           MUST BE PRESENT IF INPUT FILE VG IS PRESENT (IFLD=2).               D3                                                  
C                                                                                                                                   
CINPUT PARAMETERS...                                                                                                                
C                                                                               D5                                                  
C11-15 NLG   = NUMBER OF LONGITUDES IN THE OUTPUT GRID EXCLUDING THE CYCLIC     D5                                                  
C              LONGITUDE. IGNORED FOR INPUT FIELDS OF KIND 'ZONL'.              D5                                                  
C16-20 NLAT  = NUMBER OF LATITUDES IN THE OUTPUT GRID.                          D5                                                  
C21-25 KIND  = OUTPUT GRID AND INTERPOLATION TYPE:                              D5                                                  
C            = 1 BILINEAR INTERPOLATION TO UNSHIFTED LAT-LON GRID THAT          D5                                                  
C                INCLUDES THE POLES AND CYCLIC LONGITUDE.                       D5                                                  
C            =-1 BILINEAR INTERPOLATION TO SHIFTED LAT-LON GRID.                D5                                                  
C            = 2 BILINEAR INTERPOLATION TO GAUSSIAN GRID THAT                   D5                                                  
C                INCLUDES CYCLIC LONGITUDE.                                     D5                                                  
C            = 3 CUBIC  INTERPOLATION TO UNSHIFTED LAT-LON GRID THAT            D5                                                  
C                INCLUDES THE POLES AND CYCLIC LONGITUDE.                       D5                                                  
C            =-3 CUBIC  INTERPOLATION TO SHIFTED LAT-LON GRID.                  D5                                                  
C            = 4 CUBIC  INTERPOLATION TO GAUSSIAN GRID THAT                     D5                                                  
C                INCLUDES CYCLIC LONGITUDE.                                     D5                                                  
C            = 5 AREA-WEIGHTED INTERPOLATION TO UNSHIFTED LAT-LON GRID THAT     D5                                                  
C                INCLUDES THE POLES AND CYCLIC LONGITUDE.                       D5                                                  
C            =-5 AREA-WEIGHTED INTERPOLATION TO SHIFTED LAT-LON GRID.           D5                                                  
C            = 6 AREA WEIGHTED INTERPOLATION TO GAUSSIAN GRID THAT              D5                                                  
C                INCLUDES CYCLIC LONGITUDE.                                     D5                                                  
C                NOTE:                                                          D5                                                  
C                LONG-UNSHIFTED GRID (SHFTLG=0) INCLUDES CYCLIC LONGITUDE.      D5                                                  
C                LONG-SHIFTED (SHFTLG!=0) GRID DOESN'T INCLUDE CYCLIC LONGITUDE.D5                                                  
C            = 0 DEFAULTS TO 1.                                                 D5                                                  
C26-30 NPKGG = OUTPUT PACKING DENSITY.                                          D5                                                  
C            = 0 USE PACKING DENSITY OF INPUT FIELDS, EXCEPT WHEN THEY PACKED   D5                                                  
C                WITH PACKING DENSITY 0 IN WHICH CASE NPKGG IS SET TO 1.        D5                                                  
C31-40 SHFTLT=   GRID SHIFT ALONG LATITUDE FROM THE POLES IN DEGREES.           D5                                                  
C41-50 SHFTLG=   GRID SHIFT ALONG LONGITUDE FROM THE GREENWICH IN DEGREES.      D5                                                  
C51-55 ISLBL = 1 PROPAGATE SUPERLABELS FROM INPUT FILE INTO OUTPUT FILE.        D5                                                  
C                OTHERWISE DO NOT OUTPUT THEM.                                  D5                                                  
C56-60 IFLD  = INPUT FIELD TYPE:                                                D5                                                  
C            = 0 INPUT FILE CONTAINS SCALAR FIELDS ON A GAUSSIAN GRID           D5                                                  
C            = 1 INPUT FILE CONTAINS U- OR V-VECTOR COMPONENTS.                 D5                                                  
C            = 2 TWO INPUT FILES CONTAIN U- AND V-VECTOR COMPONENTS.            D5                                                  
C                THE RESULTS WILL BE SAVED IN TWO OUTPUT FILES.                 D5                                                  
C            =   IFLD+10 INPUT FILE IS ON AN UNSHIFTED LAT/LON GRID.            D5                                                  
C                        POLES ARE INCLUDED. THERE IS A CYCLIC MERIDIAN.        D5                                                  
C            =   IFLD+20 INPUT FILE IS ON A LAT-SHIFTED LAT/LON GRID.           D5                                                  
C                        POLES ARE EXCLUDED. THERE IS A CYCLIC MERIDIAN.        D5                                                  
C            =   IFLD+30 INPUT FILE IS ON A LON-SHIFTED LAT/LON GRID.           D5                                                  
C                        POLES ARE INCLUDED. THERE IS NO CYCLIC MERIDIAN.       D5                                                  
C            =   IFLD+40 INPUT FILE IS ON A SHIFTED LAT/LON GRID.               D5                                                  
C                        POLES ARE EXCLUDED. THERE IS NO CYCLIC MERIDIAN.       D5                                                  
C                                                                               D5                                                  
C      NOTE:   INPUT PARAMETERS SHFTLT,SHFTLG ARE USED ONLY WHEN KIND<0.        D5                                                  
C                                                                               D5                                                  
CEXAMPLES OF INPUT CARD...                                                      D5                                                  
C                                                                               D5                                                  
C* GGALL     72   37    1    2                                                  D5                                                  
C(LINEAR INTERPOLATION OF SCALAR FIELDS FROM GAUSSIAN GRID TO UNSHIFTED         D5                                                  
C 73X37 LAT-LON GRID. OUTPUT GRIDS INCLUDE POLES AND CYCLIC LONGITUDE.          D5                                                  
C OUTPUT PACKING DENSITY 2. SKIP SUPERLABELS.)                                  D5                                                  
C                                                                               D5                                                  
C* GGALL    144   73    1                             1                         D5                                                  
C(LINEAR INTERPOLATION OF SCALAR FIELDS TO UNSHIFTED LAT-LON GRID.              D5                                                  
C USE PACKING DENSITY OF INPUT FIELDS. PROPAGATE SUPERLABELS.                   D5                                                  
C OUTPUT LAT-LON GRID SPANS THE POLES AND HAS CYCLIC LONGITUDE.)                D5                                                  
C                                                                               D5                                                  
C* GGALL    144   73    1                                  1                    D5                                                  
C(LINEAR INTERPOLATION OF A VECTOR COMPONENT TO UNSHIFTED LAT-LON GRID.         D5                                                  
C USE PACKING DENSITY OF INPUT FIELDS. SKIP SUPERLABELS.                        D5                                                  
C OUTPUT LAT-LON GRID SPANS THE POLES AND HAS CYCLIC LONGITUDE.)                D5                                                  
C                                                                               D5                                                  
C* GGALL    144   73    1                                  2                    D5                                                  
C(LINEAR INTERPOLATION OF (U,V)-VECTOR TO UNSHIFTED LAT-LON GRID.               D5                                                  
C USE PACKING DENSITY OF INPUT FIELDS. SKIP SUPERLABELS.                        D5                                                  
C OUTPUT LAT-LON GRID SPANS THE POLES AND HAS CYCLIC LONGITUDE.)                D5                                                  
C                                                                               D5                                                  
C* GGALL    144   73    5                                                       D5                                                  
C(AREA-WEIGHTED INTERPOLATION OF SCALAR FIELDS TO UNSHIFTED LAT-LON GRID.       D5                                                  
C USE PACKING DENSITY OF INPUT FIELDS.                                          D5                                                  
C OUTPUT LAT-LON GRID SPANS THE POLES AND HAS CYCLIC LONGITUDE.)                D5                                                  
C                                                                               D5                                                  
C* GGALL    128   64    2                                                       D5                                                  
C(LINEAR INTERPOLATION OF SCALAR FIELDS TO OUTPUT GAUSSIAN GRID.                D5                                                  
C USE PACKING DENSITY OF INPUT FIELDS. SKIP SUPERLABELS.)                       D5                                                  
C                                                                               D5                                                  
C* GGALL     72   36   -1    1       2.5       2.5                              D5                                                  
C(LINEAR INTERPOLATION OF SCALAR FIELDS TO SHIFTED LAT-LON GRID.                D5                                                  
C OUTPUT PACKING DENSITY 1. SKIP SUPERLABELS.                                   D5                                                  
C THE OUTPUT LAT-LON GRID DOES NOT HAVE CYCLIC LONGITUDE.)                      D5                                                  
C                                                                               D5                                                  
C* GGALL    144   72   -1    1      1.25      1.25         2                    D5                                                  
C(LINEAR INTERPOLATION OF (U,V) VECTORS TO SHIFTED LAT-LON GRID.                D5                                                  
C USE PACKING DENSITY OF INPUT FIELDS. SKIP SUPERLABELS.                        D5                                                  
C THE OUTPUT LAT-LON GRID DOES NOT HAVE A CYCLIC LONGITUDE.)                    D5                                                  
C                                                                               D5                                                  
C* GGALL    180   90   -1    2        1.        1.                              D5                                                  
C(LINEAR INTERPOLATION A SCALAR FIELD TO SHIFTED LAT-LON GRID.                  D5                                                  
C PACKING DENSITY IS 2. SKIP SUPERLABELS.                                       D5                                                  
C THE OUTPUT LAT-LON GRID DOES NOT HAVE CYCLIC LONGITUDE)                       D5                                                  
C                                                                               D5                                                  
C* GGALL    180   90   -5    2        1.        0.                              D5                                                  
C(AREA-WEIGHTED INTERPOLATION TO SHIFTED LAT-LON GRID.                          D5                                                  
C OUTPUT GRID IS SHIFTED BY 1 DEGREES FROM THE POLES AND HAS CYCLIC LONGITUDE.) D5                                                  
C                                                                               D5                                                  
C* GGALL    128   64    6                             1                         D5                                                  
C(AREA-WEIGHTED INTERPOLATION TO 128X64 GAUSSIAN GRID.                          D5                                                  
C PROPAGATE SUPERLABELS.)                                                       D5                                                  
C                                                                               D5                                                  
C* GGALL     72   36   -5    0       2.5       2.5    1   40                    D5                                                  
C(AREA-WEIGHTED INTERPOLATION TO 72x36 SHIFTED LAT/LON GRID.                    D5                                                  
C PROPAGATE SUPERLABELS. INPUT GRID IS SHIFTED LAT/LON GRID.                    D5                                                  
C                                                                               D5                                                  
C* GGALL    360  180   -1    1       0.5       0.5    1                         D5                                                  
C(BILINEAR INTERPOLATION OF SCALAR FIELDS TO SHIFTED LAT-LON GRID.              D5                                                  
C OUTPUT PACKING DENSITY 1. PROPAGATE SUPERLABELS.                              D5                                                  
C THE OUTPUT LAT-LON GRID DOES NOT HAVE CYCLIC LONGITUDE.)                      D5                                                  
C----------------------------------------------------------------------------                                                       
C                                                                                                                                   
C      IMPLICIT REAL (A-H,O-Z),                                                                                                     
C     +INTEGER (I-N)                                                                                                                
      IMPLICIT NONE                                                                                                                 
      INTEGER NLTX,NLGX,NWDSX,NWDSX2                                                                                                
      PARAMETER(NLTX=181,NLGX=361,NWDSX=NLTX*NLGX,NWDSX2=1*NWDSX)                                                                   
      REAL GG(NWDSX),VG(NWDSX),GL(NWDSX),VL(NWDSX),FL(NWDSX)                                                                        
      LOGICAL OK,LWGHT                                                                                                              
      REAL SL(NLTX),CL(NLTX),WL(NLTX),WOSSL(NLTX),RAD(NLTX)                                                                         
      REAL SLAT (NLTX),SLON (NLGX),DLAT (NLTX),DLON (NLGX),                                                                         
     1     SLATH(NLTX),SLONH(NLGX),DLATH(NLTX),DLONH(NLGX)                                                                          
      INTEGER I,J,ILG,ILG1,ILAT,ILAT1,ILG1OLD,ILATOLD,                                                                              
     1     NLG,NLG1,NLAT,NLAT1,IJ,IJ1,                                                                                              
     2     KIND,NPKGG,ISLBL,INTERP,IFLD,IFLDR,IFLDM,ICYC,                                                                           
     3     NR,NFF,NOUT,NOUTW,LLG                                                                                                    
      REAL DX,DY,DXH,DYH,SHFTLG,SHFTLT,R2D,D2R,WLT,ASIN                                                                             
C                                                                                                                                   
      INTEGER IBUF,IDAT,MAXX,NC4TO8                                                                                                 
      COMMON/ICOM/IBUF(8),IDAT(NWDSX2)                                                                                              
      INTEGER JBUF(8)                                                                                                               
      DATA MAXX/NWDSX2/                                                                                                             
      REAL PI,SPVAL,SPVALT,SMALL                                                                                                    
      DATA PI/3.14159265358979E0/,SPVAL/1.E+38/,SPVALT/1.E+32/,                                                                     
     1     SMALL/1.E-6/                                                                                                             
C----------------------------------------------------------------------                                                             
                                                                                                                                    
      R2D=180.E0/PI                                                                                                                 
      D2R=PI/180.E0                                                                                                                 
                                                                                                                                    
      NFF=7                                                                                                                         
      CALL JCLPNT(NFF,11,12,13,14,15,5,6)                                                                                           
      NFF=NFF-2                                                                                                                     
      REWIND 11                                                                                                                     
      REWIND 12                                                                                                                     
C                                                                                                                                   
C     * READ INPUT CARD                                                                                                             
C                                                                                                                                   
      READ(5,5010,END=900) NLG,NLAT,KIND,NPKGG,SHFTLT,SHFTLG,                   D4                                                  
     1                     ISLBL,IFLD                                           D4                                                  
C                                                                                                                                   
C     * CHECK INTERPOLATION TYPE: 1=BILINEAR (DEFAULT), 3=CUBIC, 5=AREA-WEIGHTED                                                    
C                                                                                                                                   
      IF (KIND.EQ.0)KIND=1                                                                                                          
      IF (IABS(KIND).EQ.3.OR.KIND.EQ.4) THEN                                                                                        
        INTERP=3                                                                                                                    
        WRITE(6,'(A)')'CUBIC INTERPOLATION.'                                                                                        
      ELSE IF (IABS(KIND).EQ.5.OR.KIND.EQ.6) THEN                                                                                   
        INTERP=5                                                                                                                    
        WRITE(6,'(A)')'AREA-WEIGHTED INTERPOLATION.'                                                                                
      ELSE IF (IABS(KIND).EQ.1.OR.KIND.EQ.2) THEN                                                                                   
        INTERP=1                                                                                                                    
        WRITE(6,'(A)')'BILINEAR INTERPOLATION.'                                                                                     
      ELSE                                                                                                                          
        WRITE(6,'(A)')'***ERROR: INVALID INTERPOLATION KIND.'                                                                       
        CALL                                       XIT('GGALL',-1)                                                                  
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * PRINT OUT THE CARD                                                                                                          
C                                                                                                                                   
      WRITE(6,6010) NLG,NLAT,INTERP,NPKGG,SPVAL,                                                                                    
     1     SHFTLT,SHFTLG,ISLBL,IFLD                                                                                                 
C                                                                                                                                   
C     * CHECK OUTPUT GRID DIMENSIONS                                                                                                
C                                                                                                                                   
      IF(NLG.GT.NLGX.OR.ABS(NLAT).GT.NLTX)THEN                                                                                      
        WRITE(6,*)'***ERROR: OUTPUT DIMENSIONS ARE TOO LARGE.'                                                                      
        WRITE(6,*)'NLG=',NLG,' NLGX=',NLGX                                                                                          
        WRITE(6,*)'NLAT=',NLAT,' NLTX=',NLTX                                                                                        
        CALL                                       XIT('GGALL',-2)                                                                  
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * CHECK OUTPUT GRID TYPE                                                                                                      
C                                                                                                                                   
      IF(KIND.EQ.2.OR.KIND.EQ.4.OR.KIND.EQ.6)THEN                                                                                   
C                                                                                                                                   
C       * GAUSSIAN OUTPUT GRIDS                                                                                                     
C                                                                                                                                   
        WRITE(6,'(A)')'OUTPUT GRID IS GAUSSIAN.'                                                                                    
        LLG=2                                                                                                                       
        ICYC=1                                                                                                                      
        SHFTLG=0.E0                                                                                                                 
      ELSE IF(KIND.EQ.1.OR.KIND.EQ.3.OR.KIND.EQ.5)THEN                                                                              
C                                                                                                                                   
C       * UNSHIFTED LAT/LON OUTPUT GRIDS                                                                                            
C                                                                                                                                   
        WRITE(6,'(A)')'OUTPUT GRID IS UNSHIFTED LAT/LON.'                                                                           
        LLG=1                                                                                                                       
        ICYC=1                                                                                                                      
        SHFTLT=0.E0                                                                                                                 
        SHFTLG=0.E0                                                                                                                 
      ELSE                                                                                                                          
C                                                                                                                                   
C       * SHIFTED LAT/LON OUTPUT GRIDS                                                                                              
C                                                                                                                                   
        WRITE(6,'(A)')'OUTPUT GRID IS SHIFTED LAT/LON.'                                                                             
        LLG=1                                                                                                                       
        IF(SHFTLG.EQ.0.E0) THEN                                                                                                     
          ICYC=1                                                                                                                    
        ELSE                                                                                                                        
          ICYC=0                                                                                                                    
        ENDIF                                                                                                                       
      ENDIF                                                                                                                         
      NLG1=NLG+ICYC                                                                                                                 
C                                                                                                                                   
C     * CYCLICITY IN OUPUT GRID                                                                                                     
C                                                                                                                                   
      IF(ICYC.EQ.1) THEN                                                                                                            
        WRITE(6,'(A)')                                                                                                              
     1       'OUTPUT GRID IS ASSUMED TO HAVE A CYCLIC LONGITUDE.'                                                                   
      ELSE                                                                                                                          
        WRITE(6,'(A)')                                                                                                              
     1       'OUTPUT GRID IS ASSUMED TO HAVE NO CYCLIC LONGITUDE.'                                                                  
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * SUPERLABEL TREATMENT                                                                                                        
C                                                                                                                                   
      IF(ISLBL.EQ.1)THEN                                                                                                            
        WRITE(6,'(A)')                                                                                                              
     1       'SUPERLABELS FROM INPUT FILE WILL BE WRITTEN OUT.'                                                                     
      ELSE                                                                                                                          
        WRITE(6,'(A)')                                                                                                              
     1       'SUPERLABELS FROM INPUT FILE WON''T BE WRITTEN OUT.'                                                                   
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * INPUT FIELD TYPE PARAMETER                                                                                                  
C                                                                                                                                   
      IFLDR=MOD(IFLD,10)                                                                                                            
      IFLDM=IFLD-IFLDR                                                                                                              
      IF(IFLDR.EQ.2) THEN                                                                                                           
        WRITE(6,'(A)')'INTERPOLATE A VECTOR FIELD (U,V).'                                                                           
      ELSE IF(IFLDR.EQ.1) THEN                                                                                                      
        WRITE(6,'(A)')'INTERPOLATE A VECTOR COMPONENT U OR V.'                                                                      
      ELSE IF(IFLDR.EQ.0) THEN                                                                                                      
        WRITE(6,'(A)')'INTERPOLATE A SCALAR FIELD.'                                                                                 
      ELSE                                                                                                                          
        WRITE(6,'(A)')'***ERROR: INVALID INPUT FIELD TYPE IFLD.'                                                                    
        CALL                                       XIT('GGALL',-5)                                                                  
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * CHECK IF IFLD IS INPUT FIELD PARAMETER IS CONSISTENT WITH                                                                   
C     * THE NUMBER OF ARGUMENTS                                                                                                     
C                                                                                                                                   
      LWGHT=.FALSE.                                                                                                                 
      IF(IFLDR.EQ.2)THEN                                                                                                            
        NOUT=13                                                                                                                     
        IF(NFF.EQ.4)THEN                                                                                                            
          REWIND 13                                                                                                                 
          REWIND 14                                                                                                                 
        ELSE IF(NFF.EQ.5)THEN                                                                                                       
          WRITE(6,'(A)')                                                                                                            
     1         'SAVE FILE WITH AREA FRACTION OF NON-MISSING VALUES.'                                                                
          LWGHT=.TRUE.                                                                                                              
          NOUTW=15                                                                                                                  
          REWIND 15                                                                                                                 
        ELSE                                                                                                                        
          WRITE(6,'(A)')                                                                                                            
     1         '***ERROR: WRONG NUMBER OF ARGUMENTS FOR IFLD=2.'                                                                    
          CALL                                     XIT('GGALL',-6)                                                                  
        ENDIF                                                                                                                       
      ELSE                                                                                                                          
        NOUT=12                                                                                                                     
        IF(NFF.EQ.2)THEN                                                                                                            
          REWIND 12                                                                                                                 
        ELSE IF(NFF.EQ.3)THEN                                                                                                       
          WRITE(6,'(A)')                                                                                                            
     1         'SAVE FILE WITH AREA FRACTION OF NON-MISSING VALUES.'                                                                
          LWGHT=.TRUE.                                                                                                              
          NOUTW=13                                                                                                                  
          REWIND 13                                                                                                                 
        ELSE                                                                                                                        
          WRITE(6,'(A)')                                                                                                            
     1         '***ERROR: WRONG NUMBER OF ARGUMENTS FOR IFLD=1.'                                                                    
          CALL                                     XIT('GGALL',-7)                                                                  
        ENDIF                                                                                                                       
      ENDIF                                                                                                                         
      IF(LWGHT.AND.INTERP.NE.5)THEN                                                                                                 
        WRITE(6,'(2A)')'***ERROR: OUTPUT FILE WITH AREA FRACTIONS',                                                                 
     1       ' CAN BE USED ONLY WITH AREA-WEIGHTED INTERPOLATION.'                                                                  
        CALL                                       XIT('GGALL',-8)                                                                  
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * COMPUTE OUTPUT LATITUDES IN DEGREES FROM THE EQUATOR                                                                        
C                                                                                                                                   
      IF (LLG.EQ.1)THEN                                                                                                             
        DY=(180.E0-2.E0*SHFTLT)/FLOAT(NLAT-1)                                                                                       
        DO J=1,(NLAT+1)/2                                                                                                           
          DLAT(J)=-90.E0+FLOAT(J-1)*DY+SHFTLT                                                                                       
          DLAT(NLAT-J+1)=-DLAT(J)                                                                                                   
        ENDDO                                                                                                                       
      ELSE                                                                                                                          
        CALL GAUSSG(NLAT/2,SL,WL,CL,RAD,WOSSL)                                                                                      
        CALL TRIGL (NLAT/2,SL,WL,CL,RAD,WOSSL)                                                                                      
        DO J=1,NLAT                                                                                                                 
          DLAT(J)=RAD(J)*R2D                                                                                                        
        ENDDO                                                                                                                       
        DY=DLAT(NLAT/2+1)-DLAT(NLAT/2)                                                                                              
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * COMPUTE OUTPUT LONGITUDES IN DEGREES FROM THE GREENWICH.                                                                    
C                                                                                                                                   
      DX=360.E0/FLOAT(NLG)                                                                                                          
      DO I=1,NLG1                                                                                                                   
        DLON(I)=FLOAT(I-1)*DX+SHFTLG                                                                                                
      ENDDO                                                                                                                         
      IF(INTERP.EQ.5)THEN                                                                                                           
        IF(LLG.EQ.1)THEN                                                                                                            
C                                                                                                                                   
C         * HALF LATITUDES OF OUTPUT LAT/LON GRID                                                                                   
C                                                                                                                                   
          DYH=DY/2.E0                                                                                                               
          DO J=1,NLAT                                                                                                               
            DLATH(J)=MAX(-90.E0,DLAT(J)-DYH)                                                                                        
          ENDDO                                                                                                                     
          DLATH(NLAT+1)=MIN(90.E0,DLAT(NLAT)+DYH)                                                                                   
        ELSE                                                                                                                        
C                                                                                                                                   
C         * HALF LATITUDES OF OUTPUT GAUSSIAN GRID                                                                                  
C                                                                                                                                   
          NLAT1=NLAT+1                                                                                                              
          DLATH(1)=-90.E0                                                                                                           
          DLATH(NLAT/2+1)=0.E0                                                                                                      
          WLT=-1.E0                                                                                                                 
          DO J=2,NLAT/2                                                                                                             
            WLT=WLT+WL(J-1)                                                                                                         
            DLATH(J)=ASIN(WLT)*R2D                                                                                                  
            DLATH(NLAT1-J+1)=-DLATH(J)                                                                                              
          ENDDO                                                                                                                     
          DLATH(NLAT+1)=90.E0                                                                                                       
        ENDIF                                                                                                                       
C                                                                                                                                   
C       * HALF LONGITUDES OF OUTPUT GRID                                                                                            
C                                                                                                                                   
        DXH=DX/2.E0                                                                                                                 
        DO I=1,NLG1                                                                                                                 
          DLONH(I)=DLON(I)-DXH                                                                                                      
        ENDDO                                                                                                                       
        DLONH(NLG1+1)=DLON(NLG1)+DXH                                                                                                
      ENDIF                                                                                                                         
      WRITE(6,'(A,F10.5)')'DX=',DX                                                                                                  
      WRITE(6,'(A/1000(10F10.5/))')'OUTPUT LON=',                                                                                   
     1     (DLON(I),I=1,NLG1)                                                                                                       
      WRITE(6,'(A,F10.5)')'DY=',DY                                                                                                  
      WRITE(6,'(A/1000(10F10.5/))')'OUTPUT LAT=',                                                                                   
     1     (DLAT(J),J=1,NLAT)                                                                                                       
C                                                                                                                                   
C     * READ NEXT GRID/ZONL                                                                                                         
C                                                                                                                                   
      ILG1OLD=-1                                                                                                                    
      ILATOLD=-1                                                                                                                    
      NR=0                                                                                                                          
 100  CALL GETFLD2(11,GG,-1,-1,-1,-1,IBUF,MAXX,OK)                                                                                  
      IF(.NOT.OK)THEN                                                                                                               
        IF(NR.EQ.0) CALL                           XIT('GGALL',-9)                                                                  
        CALL PRTLAB(IBUF)                                                                                                           
        WRITE(6,6020) NR                                                                                                            
        CALL                                       XIT('GGALL',0)                                                                   
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * (OPTIONAL) READ V-COMPONENT                                                                                                 
C                                                                                                                                   
      IF(IFLDR.EQ.2)THEN                                                                                                            
        DO I=1,8                                                                                                                    
          JBUF(I)=IBUF(I)                                                                                                           
        ENDDO                                                                                                                       
        CALL GETFLD2(12,VG,-1,-1,-1,-1,IBUF,MAXX,OK)                                                                                
        IF(.NOT.OK) CALL                           XIT('GGALL',-10)                                                                 
C                                                                                                                                   
C       * MAKE SURE THAT ALL RECORDS HAVE THE SAME KIND AND DIMENSIONS                                                              
C                                                                                                                                   
        CALL CMPLBL(0,IBUF,0,JBUF,OK)                                                                                               
        IF (.NOT.OK) THEN                                                                                                           
          CALL PRTLAB (IBUF)                                                                                                        
          CALL PRTLAB (JBUF)                                                                                                        
          CALL                                     XIT('GGALL',-11)                                                                 
        ENDIF                                                                                                                       
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * CHECK INPUT GRID DIMENSIONS                                                                                                 
C                                                                                                                                   
      IF(IBUF(5).GT.NLGX.OR.IBUF(6).GT.NLTX)THEN                                                                                    
        WRITE(6,*)'***ERROR: INPUT DIMENSIONS ARE TOO LARGE.'                                                                       
        WRITE(6,*)'IBUF(5)=',IBUF(5),' NLGX=',NLGX                                                                                  
        WRITE(6,*)'IBUF(6)=',IBUF(6),' NLTX=',NLTX                                                                                  
        CALL                                       XIT('GGALL',-12)                                                                 
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * GET INPUT GRID DIMENSIONS                                                                                                   
C                                                                                                                                   
      IF(IBUF(1).EQ.NC4TO8("GRID"))THEN                                                                                             
        ILG1=IBUF(5)                                                                                                                
        ILG=ILG1-1                                                                                                                  
        ILAT=IBUF(6)                                                                                                                
        IF(IFLDM.EQ.30.OR.IFLDM.EQ.40)THEN                                                                                          
C                                                                                                                                   
C         * INCREASE ILG BY ONE FOR LON-SHIFTED GRIDS                                                                               
C                                                                                                                                   
          ILG=ILG1                                                                                                                  
          ILG1=ILG+1                                                                                                                
        ENDIF                                                                                                                       
      ELSE IF(IBUF(1).EQ.NC4TO8("ZONL"))THEN                                                                                        
C                                                                                                                                   
C       * MAKE GRID FROM ZONL                                                                                                       
C                                                                                                                                   
        ILG=1                                                                                                                       
        ILG1=2                                                                                                                      
        ILAT=IBUF(5)                                                                                                                
        DO J=ILAT,1,-1                                                                                                              
          GG(J*2  )=GG(J)                                                                                                           
          GG(J*2-1)=GG(J)                                                                                                           
        ENDDO                                                                                                                       
        IF(IFLDR.EQ.2) THEN                                                                                                         
          DO J=ILAT,1,-1                                                                                                            
            VG(J*2  )=VG(J)                                                                                                         
            VG(J*2-1)=VG(J)                                                                                                         
          ENDDO                                                                                                                     
        ENDIF                                                                                                                       
        NLG=1                                                                                                                       
        NLG1=1                                                                                                                      
      ELSE                                                                                                                          
        IF(IBUF(1).EQ.NC4TO8("LABL").AND.ISLBL.EQ.1)THEN                                                                            
C                                                                                                                                   
C         * SAVE SUPERLABELS, IF REQUESTED                                                                                          
C                                                                                                                                   
          CALL PUTFLD2(NOUT,GG,IBUF,MAXX)                                                                                           
          IF(LWGHT)THEN                                                                                                             
            CALL PUTFLD2(NOUTW,GG,IBUF,MAXX)                                                                                        
          ENDIF                                                                                                                     
          IF(IFLDR.EQ.2)THEN                                                                                                        
            CALL PUTFLD2(NOUT+1,VG,IBUF,MAXX)                                                                                       
          ENDIF                                                                                                                     
        ENDIF                                                                                                                       
C                                                                                                                                   
C       * SKIP ALL OTHER RECORD KINDS                                                                                               
C                                                                                                                                   
        GOTO 100                                                                                                                    
      ENDIF                                                                                                                         
      IF(NR.EQ.0) THEN                                                                                                              
        IF(IFLDR.EQ.2)CALL PRTLAB(JBUF)                                                                                             
        CALL PRTLAB(IBUF)                                                                                                           
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * COMPUTE LATITUDES OF INPUT GRID IN DEGREES FROM THE EQUATOR                                                                 
C                                                                                                                                   
      IF (ILAT.NE.ILATOLD) THEN                                                                                                     
        IF(IFLDM.EQ.0)THEN                                                                                                          
          WRITE(6,*)'INPUT GRID IS GAUSSIAN: ILAT=',ILAT                                                                            
          IF(MOD(ILAT,2).NE.0)THEN                                                                                                  
            WRITE(6,*)'*** GRID IS GAUSSIAN BUT ILAT IS ODD.'                                                                       
            CALL                                   XIT('GGALL',-13)                                                                 
          ENDIF                                                                                                                     
C                                                                                                                                   
C         * GAUSSIAN INPUT LATITUDES                                                                                                
C                                                                                                                                   
          CALL GAUSSG(ILAT/2,SL,WL,CL,RAD,WOSSL)                                                                                    
          CALL TRIGL (ILAT/2,SL,WL,CL,RAD,WOSSL)                                                                                    
          DO J=1,ILAT                                                                                                               
            SLAT(J)=RAD(J)*R2D                                                                                                      
          ENDDO                                                                                                                     
C                                                                                                                                   
C         * HALF LATITUDES OF INPUT GAUSSIAN GRID                                                                                   
C                                                                                                                                   
          IF(INTERP.EQ.5)THEN                                                                                                       
            ILAT1=ILAT+1                                                                                                            
            SLATH(1)=-90.E0                                                                                                         
            SLATH(ILAT/2+1)=0.E0                                                                                                    
            WLT=-1.E0                                                                                                               
            DO J=2,ILAT/2                                                                                                           
              WLT=WLT+WL(J-1)                                                                                                       
              SLATH(J)=ASIN(WLT)*R2D                                                                                                
              SLATH(ILAT1-J+1)=-SLATH(J)                                                                                            
            ENDDO                                                                                                                   
            SLATH(ILAT+1)=90.E0                                                                                                     
          ENDIF                                                                                                                     
        ELSE IF(IFLDM.EQ.10.OR.IFLDM.EQ.30)THEN                                                                                     
C                                                                                                                                   
C         * LAT-UNSHIFTED LAT/LON GRIDS (POLES ARE INCLUDED)                                                                        
C                                                                                                                                   
          WRITE(6,*)'INPUT GRID IS LAT-UNSHIFTED LAT/LON: ILAT=',ILAT                                                               
          IF(MOD(ILAT,2).EQ.0)THEN                                                                                                  
            WRITE(6,*)'***WARNING: GRID IS LAT-UNSHIFTED BUT ILAT=EVEN.'                                                            
C            CALL                                  XIT('GGALL',-14)                                                                 
          ENDIF                                                                                                                     
          DY=180.E0/FLOAT(ILAT-1)                                                                                                   
          DO J=1,ILAT                                                                                                               
            SLAT(J)=-90.E0+DY*FLOAT(J-1)                                                                                            
          ENDDO                                                                                                                     
C                                                                                                                                   
C         * HALF LATITUDES OF INPUT GRID                                                                                            
C                                                                                                                                   
          IF(INTERP.EQ.5)THEN                                                                                                       
            ILAT1=ILAT+1                                                                                                            
            SLATH(1)=-90.E0                                                                                                         
            SLATH(ILAT/2+1)=0.E0                                                                                                    
            DYH=DY/2.E0                                                                                                             
            DO J=2,ILAT/2                                                                                                           
              SLATH(J)=SLAT(J-1)+DYH                                                                                                
              SLATH(ILAT1-J+1)=-SLATH(J)                                                                                            
            ENDDO                                                                                                                   
            SLATH(ILAT+1)=90.E0                                                                                                     
          ENDIF                                                                                                                     
        ELSE IF(IFLDM.EQ.20.OR.IFLDM.EQ.40)THEN                                                                                     
C                                                                                                                                   
C         * LAT-SHIFTED LAT/LON GRIDS (POLES ARE EXCLUDED)                                                                          
C                                                                                                                                   
          WRITE(6,*)'INPUT GRID IS LAT-SHIFTED: ILAT=',ILAT                                                                         
          IF(MOD(ILAT,2).EQ.1)THEN                                                                                                  
            WRITE(6,*)'*** GRID IS LAT-SHIFTED BUT ILAT IS ODD.'                                                                    
            CALL                                   XIT('GGALL',-15)                                                                 
          ENDIF                                                                                                                     
          DY=180.E0/FLOAT(ILAT)                                                                                                     
          DYH=DY/2.E0                                                                                                               
          DO J=1,ILAT                                                                                                               
            SLAT(J)=-90.E0+DY*FLOAT(J-1)+DYH                                                                                        
          ENDDO                                                                                                                     
C                                                                                                                                   
C         * HALF LATITUDES OF INPUT GRID                                                                                            
C                                                                                                                                   
          IF(INTERP.EQ.5)THEN                                                                                                       
            ILAT1=ILAT+1                                                                                                            
            SLATH(1)=-90.E0                                                                                                         
            SLATH(ILAT/2+1)=0.E0                                                                                                    
            DO J=2,ILAT/2                                                                                                           
              SLATH(J)=SLAT(J-1)+DYH                                                                                                
              SLATH(ILAT1-J+1)=-SLATH(J)                                                                                            
            ENDDO                                                                                                                   
            SLATH(ILAT+1)=90.E0                                                                                                     
          ENDIF                                                                                                                     
        ELSE                                                                                                                        
          WRITE(6,*)'*** UNKNOWN INPUT GRID TYPE.'                                                                                  
          CALL                                     XIT('GGALL',-16)                                                                 
        ENDIF                                                                                                                       
        ILATOLD=ILAT                                                                                                                
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * COMPUTE INPUT LONGITUDES IN DEGREES FROM THE GREENWICH.                                                                     
C                                                                                                                                   
      IF (ILG1.NE.ILG1OLD) THEN                                                                                                     
        IF(IFLDM.EQ.0.OR.IFLDM.EQ.10.OR.IFLDM.EQ.20)THEN                                                                            
C                                                                                                                                   
C         * UNSHIFTED LONGITUDES                                                                                                    
C                                                                                                                                   
          WRITE(6,*)'INPUT GRID IS LON-UNSHIFTED:  ILG=',ILG1                                                                       
          IF(MOD(ILG1,2).EQ.0)THEN                                                                                                  
            WRITE(6,*)'*** GRID IS LON-UNSHIFTED BUT ILG IS EVEN.'                                                                  
           CALL                                    XIT('GGALL',-17)                                                                 
          ENDIF                                                                                                                     
          DX=360.E0/FLOAT(ILG)                                                                                                      
          DO I=1,ILG1                                                                                                               
            SLON(I)=FLOAT(I-1)*DX                                                                                                   
          ENDDO                                                                                                                     
C                                                                                                                                   
C         * HALF LONGITUDES OF INPUT GRID                                                                                           
C                                                                                                                                   
          IF(INTERP.EQ.5)THEN                                                                                                       
            DXH=DX/2.E0                                                                                                             
            DO I=1,ILG1                                                                                                             
              SLONH(I)=SLON(I)-DXH                                                                                                  
            ENDDO                                                                                                                   
            SLONH(ILG1+1)=SLON(ILG1)+DXH                                                                                            
          ENDIF                                                                                                                     
        ELSE IF(IFLDM.EQ.30.OR.IFLDM.EQ.40)THEN                                                                                     
C                                                                                                                                   
C         * SHIFTED LONGITUDES                                                                                                      
C                                                                                                                                   
          WRITE(6,*)'INPUT GRID IS LON-SHIFTED:  ILG=',ILG                                                                          
          IF(MOD(ILG,2).EQ.1)THEN                                                                                                   
            WRITE(6,*)'*** GRID IS LON-UNSHIFTED BUT ILON IS ODD.'                                                                  
            CALL                                   XIT('GGALL',-18)                                                                 
          ENDIF                                                                                                                     
          DX=360.E0/FLOAT(ILG)                                                                                                      
          DXH=DX/2.E0                                                                                                               
          DO I=1,ILG1                                                                                                               
            SLON(I)=FLOAT(I-1)*DX+DXH                                                                                               
          ENDDO                                                                                                                     
C                                                                                                                                   
C         * HALF LONGITUDES OF INPUT GRID                                                                                           
C                                                                                                                                   
          IF(INTERP.EQ.5)THEN                                                                                                       
            DO I=1,ILG1                                                                                                             
              SLONH(I)=SLON(I)-DXH                                                                                                  
            ENDDO                                                                                                                   
            SLONH(ILG1+1)=SLON(ILG1)+DXH                                                                                            
          ENDIF                                                                                                                     
        ENDIF                                                                                                                       
        ILG1OLD=ILG1                                                                                                                
      ENDIF                                                                                                                         
      IF(NR.EQ.0)THEN                                                                                                               
        WRITE(6,'(A/1000(10F10.5/))')'INPUT LAT=',                                                                                  
     1       (SLAT(J),J=1,ILAT)                                                                                                     
        WRITE(6,'(A/1000(10F10.5/))')'INPUT LON=',                                                                                  
     1       (SLON(I),I=1,ILG1)                                                                                                     
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * ADD CYCLIC MERIDIAN                                                                                                         
C                                                                                                                                   
      IF(IFLDM.EQ.30.OR.IFLDM.EQ.40)THEN                                                                                            
        DO J=ILAT,1,-1                                                                                                              
          DO I=ILG,1,-1                                                                                                             
            IJ=(J-1)*ILG+I                                                                                                          
            IJ1=(J-1)*ILG1+I                                                                                                        
            GG(IJ1)=GG(IJ)                                                                                                          
            VG(IJ1)=VG(IJ)                                                                                                          
          ENDDO                                                                                                                     
          GG(IJ1+ILG)=GG(IJ1)                                                                                                       
          VG(IJ1+ILG)=VG(IJ1)                                                                                                       
        ENDDO                                                                                                                       
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * CONVERT FROM INPUT GLOBAL GRID TO OUTPUT GLOBAL GRID                                                                        
C                                                                                                                                   
      IF(INTERP.EQ.5) THEN                                                                                                          
        CALL GGIGGW(GG,VG,ILG1,ILG,ILAT,SLON,SLAT,SLONH,SLATH,                                                                      
     1              GL,VL,NLG1,NLG,NLAT,DLON,DLAT,DLONH,DLATH,                                                                      
     2              FL,INTERP,IFLDR)                                                                                                
      ELSE                                                                                                                          
        CALL GGIGG (GG,VG,ILG1,ILG,ILAT,SLON,SLAT,                                                                                  
     1              GL,VL,NLG1,NLG,NLAT,DLON,DLAT,                                                                                  
     2              INTERP,IFLDR)                                                                                                   
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * SAVE OUTPUT GRID/ZONL                                                                                                       
C                                                                                                                                   
      IF(IBUF(1).EQ.NC4TO8("GRID"))THEN                                                                                             
        IBUF(5)=NLG1                                                                                                                
        IBUF(6)=NLAT                                                                                                                
      ELSE                                                                                                                          
        IBUF(5)=NLAT                                                                                                                
        IBUF(6)=1                                                                                                                   
      ENDIF                                                                                                                         
      IF(NPKGG.GT.0)THEN                                                                                                            
        IBUF(8)=NPKGG                                                                                                               
      ELSE                                                                                                                          
        IF(IBUF(8).EQ.0)IBUF(8)=1                                                                                                   
      ENDIF                                                                                                                         
      IF(IFLDR.EQ.2)THEN                                                                                                            
        CALL PUTFLD2(NOUT+1,VL,IBUF,MAXX)                                                                                           
        IBUF(3)=JBUF(3)                                                                                                             
      ENDIF                                                                                                                         
      CALL PUTFLD2(NOUT,GL,IBUF,MAXX)                                                                                               
C                                                                                                                                   
C     * SAVE AREA FRACTION                                                                                                          
C                                                                                                                                   
      IF(LWGHT)THEN                                                                                                                 
C        IBUF(3)=NC4TO8("WGHT")                                                                                                     
        CALL PUTFLD2(NOUTW,FL,IBUF,MAXX)                                                                                            
      ENDIF                                                                                                                         
      NR=NR+1                                                                                                                       
      GO TO 100                                                                                                                     
C                                                                                                                                   
C     * E.O.F. ON INPUT.                                                                                                            
C                                                                                                                                   
  900 CALL                                         XIT('GGALL',-19)                                                                 
C-----------------------------------------------------------------------                                                            
 5010 FORMAT(10X,4I5,2E10.0,2I5)                                                D4                                                  
 6010 FORMAT('GRID SIZE',2I6,' INTERP=',I2,' NPKGG=',I2,                                                                            
     1       ' SPVAL=',1PE12.5,' LAT/LON SHIFTS=',0P2F10.7,                                                                         
     2       ' ISLBL=',I2,' IFLD=',I2)                                                                                              
 6020 FORMAT(I6,' RECORDS PROCESSED.')                                                                                              
      END                                                                                                                           
