      PROGRAM GGAPS                                                                                                                 
C     PROGRAM GGAPS (GG,       PS,       INPUT,       OUTPUT,           )       D2                                                  
C    1         TAPE1=GG, TAPE2=PS, TAPE5=INPUT, TAPE6=OUTPUT)                                                                       
C    1         TAPE1=GG, TAPE2=PS, TAPE5=INPUT, TAPE6=OUTPUT,                                                                       
C    2         TAPE45=PLTINFO                               )                                                                       
C     -------------------------------------------------------                   D2                                                  
C                                                                               D2                                                  
C     FEB 20/08 - F.MAJAESS (ADDED "PLTINFO" FILE CREATION)                     D2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)                                                           
C     NOV 06/95 - F.MAJAESS (REVISE DEFAULT PACKING DENSITY VALUE)                                                                  
C     JUL 13/92 - E. CHAN  (DIMENSION SELECTED VARIABLES AS REAL*8)                                                                 
C     JAN 29/92 - E. CHAN  (CONVERT HOLLERITH LITERALS TO ASCII)                                                                    
C     MAY 11/83 - R.LAPRISE.                                                                                                        
C     DEC 01/80 - J.D.HENDERSON                                                                                                     
C                                                                               D2                                                  
CGGAPS   - CONVERT GAUSSIAN GRID FILE TO POLAR STEREOGRAPHIC GRID FILE  1  1 C GD1                                                  
C                                                                               D3                                                  
CAUTHOR  - J.D.HENDERSON                                                        D3                                                  
C                                                                               D3                                                  
CPURPOSE - CONVERTS A FILE OF GLOBAL GAUSSIAN GRIDS TO A FILE OF                D3                                                  
C          HEMISPHERIC POLAR STEREOGRAPHIC GRIDS. INTERPOLATION                 D3                                                  
C          IS EITHER LINEAR OR CUBIC.                                           D3                                                  
C          NOTE - IN ADDITION TO PS AND THE OUTPUT FILE, GGAPS ALSO             D3                                                  
C                 WRITES CO-ORDINATE INFORMATION TO "PLTINFO" FILE              D3                                                  
C                 (UNIT 45) TO BE SUBSEQUENTLY USED BY GGPLOT.                  D3                                                  
C            ***  THE USER SHOULD CALL GGPLOT IMMEDIATELY AFTER GGAPS,(IF       D3                                                  
C                 PRODUCING A PLOT SO DESIRED), AND BEFORE CALLING              D3                                                  
C                 GGAPS/SUBAREA AGAIN, (EXCEPT WHEN PLOTTING MULTIPLE FIELDS    D3                                                  
C                 WITH A SINGLE CALL TO GGPLOT).                                D3                                                  
C                                                                               D3                                                  
CINPUT FILES...                                                                 D3                                                  
C                                                                               D3                                                  
C      GG = GLOBAL GAUSSIAN GRIDS.                                              D3                                                  
C                                                                               D3                                                  
COUTPUT FILE...                                                                 D3                                                  
C                                                                               D3                                                  
C      PS = HEMISPHERIC POLAR STEREOGRAPHIC GRIDS                               D3                                                  
C                                                                                                                                   
CINPUT PARAMETERS...                                                                                                                
C                                                                               D5                                                  
C      LX,LY = POLAR STEREOGRAPHIC GRID DIMENSIONS.                             D5                                                  
C              LX=0; THE FIRST SIX VALUES ON THE CARD WILL ASSUME               D5                                                  
C                    THEIR DEFAULT VALUES:                                      D5                                                  
C                                                                               D5                                                  
C                      LX=51, LY=55, IP=26, JP=28, D60=3.81E5, DGRW=-10.E0      D5                                                  
C                                                                               D5                                                  
C                    WHICH WILL LEAD TO PRODUCING THE STANDARD POLAR            D5                                                  
C                    STEREOGRAPHIC PROJECTION OVER N OR S POLE BASED            D5                                                  
C                    ON CHOICE FOR NHEM WITH THE ASSOCIATED                     D5                                                  
C                     (DLAT1,DLON1)/(DLAT2/DLON2)                               D5                                                  
C                    CARTESIAN COORDINATE IN DEGREES OF                         D5                                                  
C                        LOWER_LEFT/UPPER_RIGHT                                 D5                                                  
C                    CORNER OF THE GRID BEING, RESPECTIVELY:                    D5                                                  
C                       DLAT1=-9.4., DLON1=-122.8.,                             D5                                                  
C                       DLAT2=-9.4., DLON2=57.2,                                D5                                                  
C                    IF NHEM=1 (NH), AND                                        D5                                                  
C                       DLAT1=9.4., DLON1=122.8.,                               D5                                                  
C                       DLAT2=9.4., DLON2=-57.2,                                D5                                                  
C                    IF NHEM=2 (SH).                                            D5                                                  
C                                                                               D5                                                  
C      IP,JP = GRID COORDINATES OF THE POLE                                     D5                                                  
C      D60   = GRID LENGTH AT 60 DEGREES (METERS)                               D5                                                  
C      DGRW  = ORIENTATION IN THE GRID OF GREENWICH MERIDIAN (DEG)              D5                                                  
C      NHEM  = 1 FOR N. HEM., 2 FOR S. HEM.                                     D5                                                  
C              DEFAULT VALUE IS 1.                                              D5                                                  
C      KIND  = INTERPOLATION TYPE (1 = LINEAR, OTHERWISE CUBIC).                D5                                                  
C                                                                               D5                                                  
CEXAMPLE OF CARD READ...                                                        D5                                                  
C                                                                               D5                                                  
C*   GGAPS   51   55   26   28    3.81E5      -10.    1    1                    D5                                                  
C                                                                               D5                                                  
C THIS CARD WILL PRODUCE THE STANDARD NH POLAR STEREOGRAPHIC GRID.              D5                                                  
C-------------------------------------------------------------------------                                                          
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      COMMON/BLANCK/GPS(18528),GG(18528)                                                                                            
      LOGICAL OK, EX, STAT                                                                                                          
      REAL SLAT(96)                                                                                                                 
      REAL*8 SL(96),CL(96),WL(96),WOSSL(96),RAD(96)                                                                                 
      REAL   ANG(96)                                                                                                                
C                                                                                                                                   
      COMMON/ICOM/IBUF(8),IDAT(18528)                                                                                               
      DATA MAXX/18528/                                                                                                              
C----------------------------------------------------------------------                                                             
      NFF=4                                                                                                                         
      CALL JCLPNT(NFF,1,2,5,6)                                                                                                      
      REWIND 1                                                                                                                      
      REWIND 2                                                                                                                      
C                                                                                                                                   
C     * READ P.S. GRID SIZE FROM CARD.                                                                                              
C     * KIND=1 FOR LINEAR INTERPOLATION, OTHERWISE CUBIC.                                                                           
C     * GRID IS (LX,LY) WITH POLE AT (IP,JP).                                                                                       
C     * GRID SIZE IS D60(M) AND ORIENTATION IS DGRW (DEG).                                                                          
C     * NHEM = 1,2 FOR N,S HEMISPHERE.                                                                                              
C     * BLANK CARD DEFAULTS TO STANDARD CMC GRID AND CUBIC INTERP.                                                                  
C                                                                                                                                   
      READ(5,5010,END=902) LX,LY,IP,JP,D60,DGRW,NHEM,KIND                       D4                                                  
      INTERP=3                                                                                                                      
      IF(KIND.EQ.1) INTERP=1                                                                                                        
      WRITE(6,6005) INTERP                                                                                                          
C                                                                                                                                   
      IF(NHEM.EQ.0 .OR. NHEM.EQ.3 ) NHEM=1                                                                                          
      IF(LX.GT.0) THEN                                                                                                              
          NDFLT=0                                                                                                                   
          X=FLOAT(1-IP)                                                                                                             
          Y=FLOAT(1-JP)                                                                                                             
          CALL LLFXY(DLAT1,DLON1,X,Y,D60,DGRW,NHEM)                                                                                 
          X=FLOAT(LX-IP)                                                                                                            
          Y=FLOAT(LY-JP)                                                                                                            
          CALL LLFXY(DLAT2,DLON2,X,Y,D60,DGRW,NHEM)                                                                                 
          WRITE(6,6010)LX,LY,IP,JP,D60                                                                                              
          WRITE(6,6023) DLAT1,DLON1,DLAT2,DLON2,DGRW,NHEM                                                                           
          CALL PSCAL(LY,IP,JP,D60,LX,DGRW,DLAT1,DLON1,DLAT2,DLON2,                                                                  
     1               NHEM,NDFLT,NBADPS)                                                                                             
          IF (NBADPS.NE.0) THEN                                                                                                     
            WRITE(6,6026) NBADPS                                                                                                    
            CALL                                   XIT('GGAPS',-1)                                                                  
          ENDIF                                                                                                                     
          WRITE(6,6020)LX,LY,IP,JP,D60                                                                                              
          WRITE(6,6022)DLAT1,DLON1,DLAT2,DLON2,DGRW,NHEM                                                                            
                                                                                                                                    
      ELSE                                                                                                                          
       LX=51                                                                                                                        
       LY=55                                                                                                                        
       IP=26                                                                                                                        
       JP=28                                                                                                                        
       D60=3.81E5                                                                                                                   
       DGRW=-10.E0                                                                                                                  
       IF (NHEM.EQ.1) THEN                                                                                                          
         DLAT1=-9.4E0                                                                                                               
         DLON1=-122.8E0                                                                                                             
         DLAT2=-9.4E0                                                                                                               
         DLON2=57.2E0                                                                                                               
       ELSE                                                                                                                         
         DLAT1=9.4E0                                                                                                                
         DLON1=122.8E0                                                                                                              
         DLAT2=9.4E0                                                                                                                
         DLON2=-57.2E0                                                                                                              
       ENDIF                                                                                                                        
       NDFLT=0                                                                                                                      
       WRITE(6,6010)LX,LY,IP,JP,D60                                                                                                 
       WRITE(6,6023)DLAT1,DLON1,DLAT2,DLON2,DGRW,NHEM                                                                               
      ENDIF                                                                                                                         
C                                                                                                                                   
C                                                                                                                                   
C     * CHECK AND WRITE CO-ORDINATE INFORMATION OUT FOR GGPLOT                                                                      
C                                                                                                                                   
      NOUTYP=1                                                                                                                      
      INQUIRE(FILE='pltinfo',EXIST=EX)                                                                                              
      INQUIRE(45,OPENED=STAT)                                                                                                       
      IF (STAT) THEN                                                                                                                
       CLOSE(45)                                                                                                                    
      ENDIF                                                                                                                         
      OPEN(45,FILE='pltinfo',FORM='UNFORMATTED')                                                                                    
      REWIND 45                                                                                                                     
      IF (EX) THEN                                                                                                                  
        CALL READPLTINFO(45,DLAT1P,DLON1P,DLAT2P,DLON2P,DGRWP,                                                                      
     1                   NHEMP,NOUTYPP,IOST)                                                                                        
        IF(IOST.LT.0) GOTO 120                                                                                                      
                                                                                                                                    
        OK = (DLAT1P.EQ.DLAT1) .AND. (DLON1P  .EQ.DLON1  ) .AND.                                                                    
     1       (DLAT2P.EQ.DLAT2) .AND. (DLON2P  .EQ.DLON2  ) .AND.                                                                    
     2       (DGRWP .EQ.DGRW ) .AND. (NHEMP   .EQ.NHEM   ) .AND.                                                                    
     3                               (NOUTYPP .EQ.NOUTYP)                                                                           
        PRINT *,' P=',DLAT1P,DLON1P,DLAT2P,DLON2P,DGRWP,NHEMP                                                                       
        PRINT *,' O=',DLAT1,DLON1,DLAT2,DLON2,DGRW,NHEM                                                                             
        IF (OK) THEN                                                                                                                
          GO TO 115                                                                                                                 
        ELSE                                                                                                                        
          WRITE(6,6040) DLAT1,DLON1,DLAT2,DLON2,DGRW,NHEM,NOUTYP                                                                    
          GO TO 130                                                                                                                 
        ENDIF                                                                                                                       
  120   WRITE(6,6035)                                                                                                               
  130   CONTINUE                                                                                                                    
        CLOSE(45)                                                                                                                   
        CALL SYSTEM('\rm -f pltinfo')                                                                                               
        OPEN(45,FILE='pltinfo',FORM='UNFORMATTED')                                                                                  
        REWIND 45                                                                                                                   
      ENDIF                                                                                                                         
C                                                                                                                                   
      CALL WRITEPLTINFO(45,DLAT1,DLON1,DLAT2,                                                                                       
     1     DLON2,DGRW,NHEM,NOUTYP,IOST)                                                                                             
C                                                                                                                                   
C     * GET GAUSSIAN GRID SIZE FROM FIRST FIELD IN THE FILE.                                                                        
C                                                                                                                                   
  115 CALL RECGET(1,NC4TO8("GRID"),-1,-1,-1,IBUF,MAXX,OK)                                                                           
      IF(.NOT.OK) CALL                             XIT('GGAPS',-2)                                                                  
      REWIND 1                                                                                                                      
      WRITE(6,6025) IBUF                                                                                                            
      ILG1=IBUF(5)                                                                                                                  
      ILAT=IBUF(6)                                                                                                                  
      NPACK=MIN(2,IBUF(8))                                                                                                          
C                                                                                                                                   
C     * COMPUTE LATITUDE VALUES IN DEGREES FROM THE SOUTH POLE.                                                                     
C                                                                                                                                   
      ILATH=ILAT/2                                                                                                                  
      CALL GAUSSG(ILATH,SL,WL,CL,RAD,WOSSL)                                                                                         
      CALL TRIGL (ILATH,SL,WL,CL,RAD,WOSSL)                                                                                         
      DO 110 I=1,ILAT                                                                                                               
      ANG(I)=RAD(I)*180.E0/3.14159E0                                                                                                
  110 SLAT(I)=ANG(I)+90.E0                                                                                                          
C                                                                                                                                   
C     * READ THE NEXT GAUSSIAN GRID.                                                                                                
C                                                                                                                                   
      NR=0                                                                                                                          
  150 CALL GETFLD2(1,GG,NC4TO8("GRID"),-1,-1,-1,IBUF,MAXX,OK)                                                                       
      IF(.NOT.OK.AND.NR.EQ.0) CALL                 XIT('GGAPS',-3)                                                                  
      IF(.NOT.OK)THEN                                                                                                               
        WRITE(6,6030) NR                                                                                                            
        CALL                                       XIT('GGAPS',0)                                                                   
      ENDIF                                                                                                                         
      IF(NR.EQ.0) THEN                                                                                                              
        CALL PRTLAB (IBUF)                                                                                                          
        NPACK=MIN(NPACK,IBUF(8))                                                                                                    
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * CONVERT TO P.S. GRID.                                                                                                       
C                                                                                                                                   
      CALL GGIPS2(GPS,LX,LY,IP,JP,D60,DGRW,NHEM,                                                                                    
     1             GG,ILG1,ILAT,SLAT,INTERP)                                                                                        
C                                                                                                                                   
C     * SAVE THE GRIDS.                                                                                                             
C                                                                                                                                   
      CALL SETLAB(IBUF,NC4TO8("GRID"),-1,-1,-1,LX,LY,NHEM,NPACK)                                                                    
      CALL PUTFLD2(2,GPS,IBUF,MAXX)                                                                                                 
      IF(NR.EQ.0) WRITE(6,6025) IBUF                                                                                                
      NR=NR+1                                                                                                                       
      GO TO 150                                                                                                                     
C                                                                                                                                   
C     * E.O.F. ON INPUT.                                                                                                            
C                                                                                                                                   
  902 CALL                                         XIT('GGAPS',-4)                                                                  
C-----------------------------------------------------------------------                                                            
 5010 FORMAT(10X,4I5,2E10.0,2I5)                                                D4                                                  
 6005 FORMAT('0INTERP =',I6)                                                                                                        
 6010 FORMAT('0LX,LY,IP,JP,D60 =',4I6,E12.4)                                                                                        
 6020 FORMAT('0RESULTING LX,LY,IP,JP,D60 =',4I5,E12.4)                                                                              
 6022 FORMAT('0ADJUSTED DLAT1,DLON1,DLAT2,DLON2,DGRW,NHEM = ',5E12.4,I5)                                                            
 6023 FORMAT('0CHOSEN DLAT1,DLON1,DLAT2,DLON2,DGRW,NHEM = ',5E12.4,I5)                                                              
 6025 FORMAT(1X,A4,I10,2X,A4,I10,4I6)                                                                                               
 6026 FORMAT(1X,'GGAPS-ERROR : SUB-AREA CO-ORDINATES INVALID. CODE ',I5)                                                            
 6030 FORMAT('0',I6,' RECORDS READ')                                                                                                
 6035 FORMAT(/' ** PLTINFO "FT45" FILE EXISTS AND EMPTY !'/)                                                                        
 6040 FORMAT(/,                                                                                                                     
     1       ' **  WARNING   ** SUB-AREA COORDINATES CHANGED IN',                                                                   
     2       ' PLTINFO "FT45" FILE FOR GGPLOT.'/                                                                                    
     3       ' **  REMINDER  ** NORMAL PROCEDURE IS TO PLOT',                                                                       
     4       ' PREVIOUSLY "SUBAREA/GGAPS/COFAPS" GENERATED FIELDS BY',                                                              
     5       ' GGPLOT BEFORE '/                                                                                                     
     6       '                  PROCEEDING WITH ANOTHER SET OF',                                                                    
     7       ' FIELDS HAVING DIFFERENT SUB-AREA COORDINATES.'//                                                                     
     8       ' ** NEW VALUES ** DLAT1,DLON1,DLAT2,DLON2,DGRW,NHEM,',                                                                
     9       'NOUTYP=',5E12.4,2I5/)                                                                                                 
      END                                                                                                                           
