      PROGRAM EXTRAC                                                                                                                
C     PROGRAM EXTRAC (XIN,       XOUT,       INPUT,       OUTPUT,       )       B2                                                  
C    1          TAPE1=XIN, TAPE2=XOUT, TAPE5=INPUT, TAPE6=OUTPUT)                                                                   
C     -----------------------------------------------------------               B2                                                  
C                                                                               B2                                                  
C     MAR 17/04 - F.MAJAESS,S.KHARIN (REVISED TO PROPERLY HANDLE THE CASE OF    B2                                                  
C                                     I1<0,I2>0 LONGITUDINAL DEGREE CASE        B2                                                  
C                                     MAPPING TO THE SAME INDICES)              B2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)                                                           
C     NOV 06/95 - F.MAJAESS  (REVISE TO ALLOW KHEM=3)                                                                               
C     JUL 21/92 - E. CHAN    (REPLACE UNFORMATTED I/O WITH I/O ROUTINES)                                                            
C     JUL 13/92 - E. CHAN    (DIMENSION SELECTED VARIABLES AS REAL*8)                                                               
C     JAN 29/92 - E. CHAN    (CONVERT HOLLERITH LITERALS TO ASCII)                                                                  
C     OCT 26/88 - F.MAJAESS. (ALLOW SPECIFYING DEGREES AS WELL AS INDICES                                                           
C                             IN THE INPUT CARD)                                                                                    
C     SEP 09/85 - B.DUGAS.                                                                                                          
C                                                                               B2                                                  
CEXTRAC  - ZEROES OUT ALL OF A GRID EXCEPT WINDOW(S) DEFINED IN INPUT   1  1 C  B1                                                  
C                                                                               B3                                                  
CAUTHOR  - B. DUGAS                                                             B3                                                  
C                                                                               B3                                                  
CPURPOSE - EXTRACTS CERTAIN PARTS OF AN ARRAY IN XIN, COMPLEX OR REAL, AND      B3                                                  
C          THEN PUTS THE RESULT IN  XOUT = MASK * XIN, WHERE                    B3                                                  
C                        _                                                      B3                                                  
C                       !  1 IF (I1 <= I <= I2) AND (J1 <=J <= J2)              B3                                                  
C            MASK(I,J) =<                                                       B3                                                  
C                       !  0 EVERYWHERE ELSE.                                   B3                                                  
C                        -                                                      B3                                                  
C          PLEASE REFER TO THE INPUT PARAMETERS SECTION BELOW.                  B3                                                  
C          NOTE - IF  KIND  IS  SPEC OR FOUR  THEN  I1 AND I2  DENOTE  THE      B3                                                  
C                 BEGINNING AND END OF THE RETAINED COMPLEX VECTOR.             B3                                                  
C                 IF KIND IS GRID (GLOBAL OR HEMISPHERIC)  THEN  THE  DATA      B3                                                  
C                 IS ASSUMED TO BE ORDERED FROM SOUTH TO NORTH.                 B3                                                  
C                 ALSO, AT LEAST ONE  AND UP TO  40 SETS OF I1, I2, J1 AND      B3                                                  
C                 J2 ARE READ IN FROM INPUT.                                    B3                                                  
C                 INPUT CARDS BETWEEN 20 AND 40 MAY GET IGNORED,(DEPENDING      B3                                                  
C                 ON THE TYPE AND  RANGE  OF  VALUES  SPECIFIED). INPUT IS      B3                                                  
C                 READ UNTIL "EOF" IS REACHED.                                  B3                                                  
C                                                                               B3                                                  
CINPUT FILE...                                                                  B3                                                  
C                                                                               B3                                                  
C      XIN  = A CCRN STANDARD FILE                                              B3                                                  
C                                                                               B3                                                  
COUTPUT FILE...                                                                 B3                                                  
C                                                                               B3                                                  
C      XOUT = MASKED VERSION OF XIN AS DESCRIBED ABOVE.                         B3                                                  
C                                                                                                                                   
CINPUT PARAMETERS...                                                                                                                
C                                                                               B5                                                  
C      I1   = VECTOR CONTAINING THE FIRST COORDINATES (OR THE LONGITUDINAL      B5                                                  
C             DEGREE, IF IFLG.NE.0) OF  THE LOWER LEFT-HAND CORNER  OF THE      B5                                                  
C             BOXES WE WISH TO RETAIN FROM THE ORIGINAL DATA.                   B5                                                  
C      I2   = VECTOR CONTAINING THE FIRST COORDINATES (OR THE LONGITUDINAL      B5                                                  
C             DEGREE, IF IFLG.NE.0) OF THE HIGHER RIGHT-HAND CORNER OF THE      B5                                                  
C             BOXES WE WISH TO RETAIN FROM THE ORIGINAL DATA.                   B5                                                  
C             (I2 >= I1).                                                       B5                                                  
C      J1   = VECTOR CONTAINING THE SECOND COORDINATES (OR THE LATITUDINAL      B5                                                  
C             DEGREE, IF IFLG.NE.0) OF  THE LOWER LEFT-HAND CORNER  OF THE      B5                                                  
C             BOXES WE WISH TO RETAIN FROM THE ORIGINAL DATA.                   B5                                                  
C      J2   = VECTOR CONTAINING THE SECOND COORDINATES (OR THE LATITUDINAL      B5                                                  
C             DEGREE, IF IFLG.NE.0) OF THE HIGHER RIGHT-HAND CORNER OF THE      B5                                                  
C             BOXES WE WISH TO RETAIN FROM THE ORIGINAL DATA.                   B5                                                  
C             (J2 >= J1).                                                       B5                                                  
C      IFLG = USED TO IDENTIFY THE VALUES SUPPLIED IN I1, I2, J1 AND J2 AS      B5                                                  
C             INDICES (IFLG.EQ.0) OR AS DEGREES (IFLG.NE.0)                     B5                                                  
C             NOTE - IFLG MUST BE ZERO FOR ALL THE FIELDS WHOSE KIND IS         B5                                                  
C                    DIFFERENT FROM GRID (I.E. I1, I2, J1 AND J2 MAY  HAVE      B5                                                  
C                    DEGREE VALUES ONLY FOR GRID FIELDS, OTHERWISE THEY         B5                                                  
C                    MUST CONTAIN INDICES VALUES)                               B5                                                  
C                    IF IFLG.NE.0 THEN                                          B5                                                  
C                       I1 , I2 VALUES (-359,360) WITH I2 >= I1.                B5                                                  
C                       J1 , J2 VALUES (- 90, 90) WITH J2 >= J1.                B5                                                  
C      NTYP = USED TO DETERMINE THE TYPE OF GRID (IGNORED  IF  IFLG.EQ.0).      B5                                                  
C             ZERO DEFAULTS TO A GAUSSIAN GRID  AND  ANYTHING  ELSE  TO  A      B5                                                  
C             LAT-LONG GRID.                                                    B5                                                  
C                                                                               B5                                                  
CEXAMPLE OF INPUT CARD(S)...                                                    B5                                                  
C                                                                               B5                                                  
C*EXTRAC     11   41    1  101                                                  B5                                                  
C*          -10   40  -90   10    1    1                                        B5                                                  
C*          271  361   31  161    0                                             B5                                                  
C*          270  360  -60   70    1    1                                        B5                                                  
C*            0   45   80   89    1    1                                        B5                                                  
C*         -330   20  -90   90    1    0                                        B5                                                  
C*            1   46  171  180                                                  B5                                                  
C--------------------------------------------------------------------------                                                         
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      REAL MASK(65341),ALON(361),ALAT(181),DLON,DLAT,                                                                               
     1     OFFLAT,API,ARADEG                                                                                                        
      REAL*8  SL(181),CL(181),WL(181),WOSSL(181),RAD(181)                                                                           
      INTEGER I1(40), I2(40), J1(40), J2(40)                                                                                        
      COMMON/BLANCK/G(65341)                                                                                                        
                                                                                                                                    
      LOGICAL OK,SPEC,DEGCNV,LATLON                                                                                                 
      COMMON/ICOM/IBUF(8),IDAT(65341)                                                                                               
      DATA MAXX/65341/                                                                                                              
C     DATA NSPEC/4HSPEC/, NFOUR/4HFOUR/ ,NGRID/4HGRID/                                                                              
C---------------------------------------------------------------------                                                              
      NSPEC=NC4TO8("SPEC")                                                                                                          
      NFOUR=NC4TO8("FOUR")                                                                                                          
      NGRID=NC4TO8("GRID")                                                                                                          
      NFF=4                                                                                                                         
      CALL JCLPNT(NFF,1,2,5,6)                                                                                                      
      REWIND 1                                                                                                                      
      REWIND 2                                                                                                                      
C                                                                                                                                   
C     * GET FIELD KIND(S) AND GRID DIMENSIONS.                                                                                      
C                                                                                                                                   
      CALL FBUFFIN(1,IBUF,-8,K,LEN)                                                                                                 
      IF (K.GE.0) GOTO 900                                                                                                          
      CALL PRTLAB(IBUF)                                                                                                             
      REWIND 1                                                                                                                      
C                                                                                                                                   
      IF((IBUF(1).EQ.NSPEC).OR.(IBUF(1).EQ.NFOUR)) THEN                                                                             
        SPEC = .TRUE.                                                                                                               
        ILG  = IBUF(5)*2                                                                                                            
      ELSE                                                                                                                          
        SPEC = .FALSE.                                                                                                              
        ILG  = IBUF(5)                                                                                                              
      ENDIF                                                                                                                         
      KIND = IBUF(1)                                                                                                                
      NLG  = IBUF(5)                                                                                                                
      NLT  = IBUF(6)                                                                                                                
      KHEM = IBUF(7)                                                                                                                
      NWDS = ILG*NLT                                                                                                                
C                                                                                                                                   
C     * SET UP FOR POSSIBLE CONVERSION FROM DEGREES TO INDICES.                                                                     
C                                                                                                                                   
      DEGCNV=((KIND.EQ.NGRID).AND.(KHEM.GE.0).AND.(KHEM.LE.3))                                                                      
C                                                                                                                                   
      IF(DEGCNV)THEN                                                                                                                
        DLON=360.E0/FLOAT(ILG-1)                                                                                                    
        DLON=1.E-10*ANINT(1.E10*DLON)                                                                                               
        DO 10 I=1,ILG                                                                                                               
          ALON(I)=FLOAT(I-1)*DLON                                                                                                   
   10   CONTINUE                                                                                                                    
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * READ I1, I2, J1, J2 AND IFLG FROM CARD.                                                                                     
C                                                                                                                                   
      NCONV=0                                                                                                                       
      N=-1                                                                      B4                                                  
      WRITE(6,6000)                                                                                                                 
  100 N=N+1                                                                     B4                                                  
      IF(N.GE.40) THEN                                                          B4                                                  
        GOTO 150                                                                B4                                                  
      ENDIF                                                                     B4                                                  
      READ(5,5010,END=150) I1(N+1),I2(N+1),J1(N+1),J2(N+1),IFLG,NTYP            B4                                                  
      WRITE(6,6005) I1(N+1),I2(N+1),J1(N+1),J2(N+1),IFLG,NTYP                                                                       
      IF ( I1(N+1).GT.I2(N+1) .OR. J1(N+1).GT.J2(N+1)) THEN                                                                         
        WRITE(6,6040)                                                                                                               
        CALL                                       XIT('EXTRAC',-1)                                                                 
      ENDIF                                                                                                                         
      IF(IFLG.NE.0)THEN                                                                                                             
        IF(.NOT.DEGCNV)CALL                        XIT('EXTRAC',-2)                                                                 
        IF(NCONV.EQ.0)THEN                                                                                                          
          IF(NTYP.NE.0)THEN                                                                                                         
            LATLON=.TRUE.                                                                                                           
            IF(KHEM.EQ.0 .OR. KHEM.EQ.3 )THEN                                                                                       
              DLAT=180.E0/FLOAT(NLT-1)                                                                                              
            ELSE                                                                                                                    
              DLAT= 90.E0/FLOAT(NLT-1)                                                                                              
            ENDIF                                                                                                                   
            DLAT=1.E-10*ANINT(1.E10*DLAT)                                                                                           
            IF(KHEM.EQ.1)THEN                                                                                                       
              OFFLAT=  0.0E0                                                                                                        
            ELSE                                                                                                                    
              OFFLAT=-90.0E0                                                                                                        
            ENDIF                                                                                                                   
            DO 120 J=1,NLT                                                                                                          
              ALAT(J)=OFFLAT+FLOAT(J-1)*DLAT                                                                                        
  120       CONTINUE                                                                                                                
          ELSE                                                                                                                      
            LATLON=.FALSE.                                                                                                          
            IF(KHEM.EQ.0 .OR. KHEM.EQ.3 )THEN                                                                                       
              NLTH=NLT/2                                                                                                            
            ELSE                                                                                                                    
              NLTH=NLT                                                                                                              
            ENDIF                                                                                                                   
            CALL GAUSSG(NLTH,SL,WL,CL,RAD,WOSSL)                                                                                    
            CALL TRIGL (NLTH,SL,WL,CL,RAD,WOSSL)                                                                                    
            API=4.0E0*ATAN(1.0E0)                                                                                                   
            ARADEG=180.E0/API                                                                                                       
            IF(KHEM.EQ.1)THEN                                                                                                       
              DO 125 J=1,NLT                                                                                                        
                ALAT(J)=RAD(NLT+J)*ARADEG                                                                                           
  125         CONTINUE                                                                                                              
            ELSE                                                                                                                    
              DO 130 J=1,NLT                                                                                                        
                ALAT(J)=RAD(J)*ARADEG                                                                                               
  130         CONTINUE                                                                                                              
            ENDIF                                                                                                                   
          ENDIF                                                                                                                     
        ENDIF                                                                                                                       
C                                                                                                                                   
        IF(((NTYP.NE.0).AND.(.NOT.LATLON)).OR.                                                                                      
     1     ((NTYP.EQ.0).AND.(     LATLON))) CALL   XIT('EXTRAC',-3)                                                                 
C                                                                                                                                   
        I1NP1=I1(N+1)                                                                                                               
        IF(I1(N+1).LT.0) I1(N+1)=I1(N+1)+360                                                                                        
        KI1=0                                                                                                                       
        DLON1=FLOAT(I1(N+1))                                                                                                        
        IF(DLON1.LE.ALON( 1 )) KI1=1                                                                                                
        IF(DLON1.GE.ALON(ILG)) KI1=ILG                                                                                              
        IF(KI1.EQ.0)THEN                                                                                                            
          KI1=MIN(ISRCHFGE(ILG,ALON,1,DLON1)-1,ILG)                                                                                 
          IF(ALON(KI1+1)-DLON1 .LT. DLON1-ALON(KI1)) KI1=KI1+1                                                                      
        ENDIF                                                                                                                       
C                                                                                                                                   
        I2NP1=I2(N+1)                                                                                                               
        IF(I2(N+1).LT.0) I2(N+1)=I2(N+1)+360                                                                                        
        KI2=0                                                                                                                       
        DLON2=FLOAT(I2(N+1))                                                                                                        
        IF(DLON2.LE.ALON( 1 )) KI2=1                                                                                                
        IF(DLON2.GE.ALON(ILG)) KI2=ILG                                                                                              
        IF(KI2.EQ.0)THEN                                                                                                            
          KI2=MIN(ISRCHFGE(ILG,ALON,1,DLON2)-1,ILG)                                                                                 
          IF(ALON(KI2+1)-DLON2 .LT. DLON2-ALON(KI2)) KI2=KI2+1                                                                      
        ENDIF                                                                                                                       
C                                                                                                                                   
        DLAT1=FLOAT(J1(N+1))                                                                                                        
        KJ1=0                                                                                                                       
        IF(DLAT1.LE.ALAT( 1 )) KJ1=1                                                                                                
        IF(DLAT1.GE.ALAT(NLT)) KJ1=NLT                                                                                              
        IF(KJ1.EQ.0)THEN                                                                                                            
          KJ1=MIN(ISRCHFGE(NLT,ALAT,1,DLAT1)-1,NLT)                                                                                 
          IF(ALAT(KJ1+1)-DLAT1 .LT. DLAT1-ALAT(KJ1)) KJ1=KJ1+1                                                                      
        ENDIF                                                                                                                       
C                                                                                                                                   
        KJ2=0                                                                                                                       
        DLAT2=FLOAT(J2(N+1))                                                                                                        
        IF(DLAT2.LT.DLAT1) KJ2=KJ1                                                                                                  
        IF(KJ2.EQ.0)THEN                                                                                                            
          IF(DLAT2.LE.ALAT( 1 )) KJ2=1                                                                                              
          IF(DLAT2.GE.ALAT(NLT)) KJ2=NLT                                                                                            
        ENDIF                                                                                                                       
        IF(KJ2.EQ.0)THEN                                                                                                            
          KJ2=MIN(ISRCHFGE(NLT,ALAT,1,DLAT2)-1,NLT)                                                                                 
          IF(ALAT(KJ2+1)-DLAT2 .LT. DLAT2-ALAT(KJ2)) KJ2=KJ2+1                                                                      
        ENDIF                                                                                                                       
C                                                                                                                                   
        IFLG=0                                                                                                                      
        IF(KI1.GT.KI2.OR.(KI1.EQ.KI2.AND.(I2NP1-I1NP1).GT.DLON))THEN                                                                
          I1(N+1)=KI1                                                                                                               
          I2(N+1)=ILG                                                                                                               
          J1(N+1)=KJ1                                                                                                               
          J2(N+1)=KJ2                                                                                                               
          WRITE(6,6005) I1(N+1),I2(N+1),J1(N+1),J2(N+1),IFLG,NTYP                                                                   
          NCONV=NCONV+1                                                                                                             
          KI1=1                                                                                                                     
          N=N+1                                                                                                                     
          IF(N.GE.40) THEN                                                                                                          
            WRITE(6,6030)                                                                                                           
            GOTO 150                                                                                                                
          ENDIF                                                                                                                     
        ENDIF                                                                                                                       
        I1(N+1)=KI1                                                                                                                 
        I2(N+1)=KI2                                                                                                                 
        J1(N+1)=KJ1                                                                                                                 
        J2(N+1)=KJ2                                                                                                                 
        WRITE(6,6005) I1(N+1),I2(N+1),J1(N+1),J2(N+1),IFLG,NTYP                                                                     
        NCONV=NCONV+1                                                                                                               
      ENDIF                                                                                                                         
      GOTO 100                                                                  B4                                                  
  150 IF (N.EQ.0) CALL                             XIT('EXTRAC',-4)                                                                 
C                                                                                                                                   
C     * ADJUST I1, I2, J1 AND J2 IF NECESSARY.                                                                                      
C                                                                                                                                   
      DO 200 L=1,N                                                                                                                  
        IF(SPEC) THEN                                                                                                               
          I1(L)=I1(L)*2-1                                                                                                           
          I2(L)=I2(L)*2                                                                                                             
        ENDIF                                                                                                                       
        I1(L)=MIN(MAX( 1   ,I1(L)), ILG )                                                                                           
        I2(L)=MAX(MIN(I2(L), ILG ),I1(L))                                                                                           
        IF(SPEC) I2(L)=MAX(I2(L),I1(L)+1)                                                                                           
        J1(L)=MIN(MAX( 1   ,J1(L)), NLT )                                                                                           
        J2(L)=MAX(MIN(J2(L), NLT ),J1(L))                                                                                           
  200 CONTINUE                                                                                                                      
C                                                                                                                                   
C     * SET MASK ARRAY.                                                                                                             
C                                                                                                                                   
      DO 300 K=1,NWDS                                                                                                               
          MASK(K)=0.0E0                                                                                                             
  300 CONTINUE                                                                                                                      
                                                                                                                                    
      DO 330 L=1,N                                                                                                                  
        DO 320 J=J1(L),J2(L)                                                                                                        
          IJ1=(J-1)*ILG + I1(L)                                                                                                     
          IJ2= IJ1 + (I2(L) - I1(L))                                                                                                
          DO 310 I=IJ1,IJ2                                                                                                          
            MASK(I)=1.0E0                                                                                                           
  310     CONTINUE                                                                                                                  
  320   CONTINUE                                                                                                                    
  330 CONTINUE                                                                                                                      
C                                                                                                                                   
C     * READ THE NEXT GRID FROM FILE XIN.                                                                                           
C                                                                                                                                   
      NR=0                                                                                                                          
  500 CALL GETFLD2(1,G,-1,-1,-1,-1,IBUF,MAXX,OK)                                                                                    
                                                                                                                                    
      IF (.NOT.OK) THEN                                                                                                             
        CALL PRTLAB(IBUF)                                                                                                           
        WRITE(6,6010) NR                                                                                                            
        CALL                                       XIT('EXTRAC',0)                                                                  
      ENDIF                                                                                                                         
                                                                                                                                    
      IF((IBUF(1).NE.KIND).OR.(IBUF(5).NE. NLG).OR.                                                                                 
     1   (IBUF(6).NE. NLT).OR.(IBUF(7).NE.KHEM)    )                                                                                
     2   CALL                                      XIT('EXTRAC',-5)                                                                 
                                                                                                                                    
      DO 600 K=1,NWDS                                                                                                               
        G(K)=G(K)*MASK(K)                                                                                                           
  600 CONTINUE                                                                                                                      
                                                                                                                                    
C     * SAVE ON FILE XOUT.                                                                                                          
                                                                                                                                    
      CALL PUTFLD2(2,G,IBUF,MAXX)                                                                                                   
      NR=NR+1                                                                                                                       
                                                                                                                                    
      GO TO 500                                                                                                                     
                                                                                                                                    
C     * PREMATURE E.O.F. ON UNIT 1.                                                                                                 
                                                                                                                                    
  900 CALL                                         XIT('EXTRAC',-6)                                                                 
C---------------------------------------------------------------------                                                              
 5010 FORMAT(10X,6I5)                                                           B4                                                  
 6000 FORMAT('0 EXTRAC',8X,'I1',8X,'I2',8X,'J1',8X,'J2',6X,'IFLG',                                                                  
     1                                                  6X,'NTYP :')                                                                
 6005 FORMAT(8X,6(1X,I9))                                                                                                           
 6010 FORMAT('0 EXTRAC READ',I6,' RECORDS.')                                                                                        
 6030 FORMAT('0 EXTRAC: NOT SUFFICIENT ARRAY SIZE TO EXPAND')                                                                       
 6040 FORMAT('0 EXTRAC: "I1 > I2" AND/OR "J1 > J2" VIOLATION!')                                                                     
      END                                                                                                                           
