      SUBROUTINE GGFILL(GLR,MASK,DLON,DLAT,LONB,LATB,ILG,ILG1,ILAT,
     1                  NBADPTS,ICHOICE,PI,SPVAL,OK)

C     * MAR 07/08 - M.LAZARE. (MODIFIED TO ADD IF/ENDIF ON NLPTS.GT.0
C     *                        WHEN ICHOICE=0, ELSE GETS MEMORY FAULT).
C     * SEP 18/06 - F.MAJAESS (MODIFIED FOR TOLERANCE IN SPVAL)
C     * MAY 28/04 - M. LAZARE. BUGFIX TO INITIALIZE NEWMAX TO ZERO BEFORE
C     *                        LOOP 150.
C     * MAY 14/04 - M. LAZARE. CHANGE THE DIMENSION OF "LONB(1),
C     *                        LATB(1)" TO "*".
C     * OCT 30/97 - M. LAZARE.

C     * SUBROUTINE TO TRY AND "BORROW" INFORMATION FROM SURROUNDING LOW-
C     * RESOLUTION GRID POINTS TO DEFINE INITIAL FIELDS AT POINTS WHERE
C     * THERE IS NO HIGH-RESOLUTION INFORMATION AVAILABLE. 
C
C     * COMPASS DIRECTIONS N,NE,E,SE,S,SW,W,NW ARE CONSIDERED (NDIR=8).
C     * ANY OF THESE POINTS WHOSE LAND MASK MATCHES THAT OF THE POINT
C     * IN QUESTION ARE USED IN AN AREA-AVERAGED SENSE. THE SUBROUTINE
C     * ABORTS IF THIS IS NOT AN AREA-AVERAGING CASE (ICHOICE=1).
C
C     * NOTE THAT IF THE VALUE OF THE SURROUNDING POINT FIELD IS "SPVAL"
C     * THIS POINT IS NOT USED IN THE CALCULATION, SINCE THIS SPECIAL
C     * VALUE IS THE ONE INITIALLY STORED IN THE FIRST PASS THROUGH
C     * THE DATA, I.E. IT HAS NOT BEEN PROCESSED YET (I.E. NORTH AND/OR
C     * EAST OF THE CURRENT PROBLEM POINT).
C
C     * IF THIS IS STILL NOT SUCCESSFUL, THE ROUTINE PASSES "OK=.FALSE."
C     * BACK UP THROUGH "HRTOLR3" TO THE PHYSICS INITIALIZATION
C     * PROGRAM, WHICH SUBSEQUENTLY ABORTS.

      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
      PARAMETER (NDIR=8)
C
C     * INPUT/OUTPUT FIELDS:
C
      REAL GLR(ILG1,ILAT), MASK(ILG1,ILAT)
      REAL DLON(ILG), DLAT(ILAT)
      INTEGER LONB(*), LATB(*)
C
C     * INTERNAL ARRAYS:
C
      REAL VAL(NDIR), DIST(NDIR)
      INTEGER NMAX(NDIR), LOCFST(NDIR), IVAL(NDIR), LOCAT(NDIR)
      INTEGER LON(NDIR), LAT(NDIR) 
C
      LOGICAL OK
C===========================================================================
C     * FIRST ABORT IF NOT AREA-AVERAGING.
C
      IF(ICHOICE.NE.0 .AND. ICHOICE.NE.1)     CALL XIT('GGFILL',-1)

C     * SETUP IN "SPVALT" THE TOLERANCE TO CHECK "SPVAL" AGAINST.

      SPVALT=1.E-6*SPVAL 
C
C     * LOOP OVER BAD POINTS.
C
      DO 300 NN=1,NBADPTS
          I=LONB(NN)
          J=LATB(NN)
C
C         * DEFINE LAT/LON INDICES OF SURROUNDING POINTS IN A CLOCKWISE
C         * SENSE (NORTH=1,NE=2,EAST=3,...NW=8).
C
          IF(J.EQ.1)            THEN
              NPTS=5   
              LAT(1)=2
              LAT(2)=2
              LAT(3)=1
              LAT(4)=1
              LAT(5)=2    
              LON(1)=I 
              IF(I.EQ.1)               THEN   
                  LON(2)=2
                  LON(3)=2   
                  LON(4)=ILG1-1
                  LON(5)=ILG1-1     
              ELSE IF(I.EQ.ILG1-1)     THEN  
                  LON(2)=1
                  LON(3)=1   
                  LON(4)=ILG1-2
                  LON(5)=ILG1-2     
              ELSE
                  LON(2)=I+1
                  LON(3)=I+1   
                  LON(4)=I-1
                  LON(5)=I-1     
              ENDIF
          ELSE IF(J.EQ.ILAT)    THEN
              NPTS=5
              LAT(1)=ILAT
              LAT(2)=ILAT-1
              LAT(3)=ILAT-1
              LAT(4)=ILAT-1
              LAT(5)=ILAT
              LON(3)=I 
              IF(I.EQ.1)               THEN   
                  LON(1)=2
                  LON(2)=2   
                  LON(4)=ILG
                  LON(5)=ILG     
              ELSE IF(I.EQ.ILG)     THEN  
                  LON(1)=1
                  LON(2)=1   
                  LON(4)=ILG-1
                  LON(5)=ILG-1     
              ELSE
                  LON(1)=I+1
                  LON(2)=I+1   
                  LON(4)=I-1
                  LON(5)=I-1     
              ENDIF
          ELSE      
              NPTS=8
              LAT(1)=J+1
              LAT(2)=J+1
              LAT(3)=J
              LAT(4)=J-1
              LAT(5)=J-1
              LAT(6)=J-1
              LAT(7)=J
              LAT(8)=J+1          
              LON(1)=I 
              LON(5)=I 
              IF(I.EQ.1)               THEN   
                  LON(2)=2
                  LON(3)=2   
                  LON(4)=2
                  LON(6)=ILG
                  LON(7)=ILG     
                  LON(8)=ILG
              ELSE IF(I.EQ.ILG)     THEN  
                  LON(2)=1
                  LON(3)=1   
                  LON(4)=1
                  LON(6)=ILG-1
                  LON(7)=ILG-1     
                  LON(8)=ILG-1 
              ELSE
                  LON(2)=I+1
                  LON(3)=I+1   
                  LON(4)=I+1
                  LON(6)=I-1     
                  LON(7)=I-1
                  LON(8)=I-1
              ENDIF     
          ENDIF
C   
          IF(ICHOICE.EQ.0 )                                   THEN
C
C 
C           * FOR DISCRETE-VALUED FIELDS, KEEP TRACK OF EACH QUALIFYING HIGH- 
C           * RESOLUTION POINT VALUE IN ARRAY VAL, AND THEIR DISTANCE FROM THE
C           * LOW-RESOLUTION GRID POINT IN ARRAY DIST.
C 
            NLPTS=0
            DO 100 ID=1,NPTS
              II=LON(ID)
              JJ=LAT(ID)
              DISTLAT=DLAT(J)-DLAT(JJ)
              FACTL=(COS(DLAT(J)*PI/180.E0))**2 
              IF(MASK(II,JJ).EQ.MASK(I,J) .AND. 
     1              ABS(GLR(II,JJ)-SPVAL).GT.SPVALT)THEN
                NLPTS=NLPTS+1
                VAL(NLPTS)=GLR(II,JJ) 
                DISTRY=ABS(DLON(I)-DLON(II))
                DISTLON=MIN(DISTRY,(360.E0-DISTRY)) 
                DIST(NLPTS)=(FACTL*(DISTLON**2)+DISTLAT**2)*
     1                      ((PI/180.E0)**2)
              ENDIF
  100       CONTINUE
C
C           * NOW DO THE REMAINING CALCULATION. 
C 
            IF(NLPTS.GT.0)                              THEN
             NEWMAX=0
             DO 150 N=1,NLPTS
              IF(N.EQ.1) THEN 
                NT=0
                NUM=0 
              ELSE
                CALL WHENEQ(N-1,VAL,1,VAL(N),IVAL,NUM)
              ENDIF 
              IF(NUM.EQ.0) THEN 
                NT=NT+1 
                CALL WHENEQ(NLPTS,VAL,1,VAL(N),IVAL,NVAL) 
                LOCFST(NT)=N
                NMAX(NT)=NVAL 
                IF(NVAL.GT.NEWMAX) NEWMAX=NVAL
              ENDIF 
  150        CONTINUE
C 
C            * IF ONLY ONE VALUE EXISTS IN THE LOW-RESOLUTION GRID SQUARE, 
C            * RETURN ITS VALUE. 
C 
             IF(NT.EQ.1) THEN
              GLR(I,J)=VAL(LOCFST(1)) 
             ELSE
C 
C             * IF ONLY ONE UNIQUE-VALUED POINT HAS THE MOST FREQUENT 
C             * OCCURRENCE (NEWMAX), RETURN ITS VALUE. 
C 
              CALL WHENEQI(NT,NMAX,1,NEWMAX,IVAL,NVAL) 
              IF(NVAL.EQ.1) THEN
                GLR(I,J)=VAL(LOCFST(IVAL(1))) 
              ELSE
C 
C               * OTHERWISE, SEARCH THROUGH THE UNIQUE-VALUED POINTS TO 
C               * DETERMINE THE VALUES HAVING THE MOST FREQUENT OCCURRENCE AS 
C               * NEWMAX AND THEIR LOCATION IN ARRAY VAL (DEFINED BY ARRAY
C               * LOCAT). 
C 
                NPP=0 
                DO 175 N=1,NT 
                  IF(NMAX(N).EQ.NEWMAX) THEN
                    TEMP=VAL(LOCFST(N)) 
                    CALL WHENEQ(NLPTS,VAL,1,TEMP,IVAL,NVAL) 
                    DO 170 IJK=1,NVAL 
                      NPP=NPP+1 
                      LOCAT(NPP)=IVAL(IJK)
  170               CONTINUE
                  ENDIF 
  175           CONTINUE
C 
C               * CHOOSE THE CLOSEST POINT TO THE TARGET HAVING THAT MAXIMUM- 
C               * OCCURRING VALUE.
C 
                DISTMIN=1.E20 
                DO 180 N=1,NPP
                  IF(DIST(LOCAT(N)).LT.DISTMIN) THEN
                    DISTMIN=DIST(LOCAT(N))
                    NCHOICE=LOCAT(N)
                  ENDIF 
  180           CONTINUE
                GLR(I,J)=VAL(NCHOICE) 
              ENDIF 
             ENDIF 
            ENDIF

          ELSE
C 
C           * FOR CONTINUOUS FIELDS, AREA-AVERAGE THESE POINTS. 
C 
            BIGAREA=0.E0
            SUM=0.E0
            NLPTS=0     
C
            DO 200 ID=1,NPTS
              II=LON(ID)
              JJ=LAT(ID)
C
              IF(JJ.EQ.ILAT) THEN 
                  SLATP=90.E0 
              ELSE
                  SLATP=0.5E0*(DLAT(JJ+1)+DLAT(JJ)) 
              ENDIF 
              IF(JJ.EQ.1) THEN
                  SLATM=-90.E0
              ELSE
                  SLATM=0.5E0*(DLAT(JJ-1)+DLAT(JJ)) 
              ENDIF 
              FACTLAT=(COS(DLAT(JJ)*PI/180.E0))*(SLATP-SLATM)*PI/180.E0 
C
              IF(II.EQ.ILG)     THEN
                  IIM1=II-1
                  SLONP=0.5E0*(DLON(II)+360.E0)
                  SLONM=0.5E0*(0.E0-DLON(ILG))
              ELSE IF(II.EQ.1)THEN
                  IIP1=II+1
                  SLONP=0.5E0*(DLON(IIP1)+DLON(II))
                  SLONM=0.5E0*(DLON(II)+0.E0)
              ELSE 
                  IIP1=II+1
                  IIM1=II-1
                  SLONP=0.5E0*(DLON(IIP1)+DLON(II))
                  SLONM=0.5E0*(DLON(IIM1)+DLON(II))
              ENDIF
C
              IF(MASK(II,JJ).EQ.MASK(I,J) .AND. 
     1              ABS(GLR(II,JJ)-SPVAL).GT.SPVALT)THEN
                  NLPTS=NLPTS+1
                  AREA=FACTLAT*((SLONP-SLONM)*PI/180.E0)
                  SUM=SUM+AREA*GLR(II,JJ) 
                  BIGAREA=BIGAREA+AREA
              ENDIF
  200       CONTINUE     
C
            IF(NLPTS.NE.0) THEN 
              GLR(I,J)=SUM/BIGAREA
            ENDIF
          ENDIF
C
C         * FINALLY, PRINT OUT RESULTS OF CALCULATION OR WARNING
C         * MESSAGE.
C
          IF(NLPTS.NE.0) THEN
              PRINT *, '0SUCCESSFUL FILL! I,J,MASK,GLR = ',I,J,
     1                   MASK(I,J),GLR(I,J)
              PRINT *, ' SURROUNDING POINTS USED WERE:'
              DO 250 ID=1,NPTS
                  II=LON(ID)
                  JJ=LAT(ID) 
                  IF(MASK(II,JJ).EQ.MASK(I,J) .AND. 
     1              ABS(GLR(II,JJ)-SPVAL).GT.SPVALT)THEN
                      PRINT *, ' I,J,MASK,GLR = ',II,JJ,
     1                           MASK(II,JJ),GLR(II,JJ)
                  ENDIF
  250         CONTINUE
          ELSE
              OK=.FALSE.
              LONBAD=I
              LATBAD=J
              PRINT *, '0UNSUCCESSFUL FILL! NO VALID SURROUNDING PTS!'
              PRINT *, ' NO POINTS FOUND WITHIN GRID SQUARE CENTRED AT',
     1                 ' I,J = ',I,J
              RETURN
          ENDIF
  300 CONTINUE
C
      RETURN
      END
