      SUBROUTINE GGIGGW(UG,VG,ILG1,ILG,ILAT,SLON,SLAT,SLONH,SLATH,
     1                  UL,VL,NLG1,NLG,NLAT,DLON,DLAT,DLONH,DLATH,
     2                  WL,INTERP,IFLD)
C
C     * MAR 25/13 - S.KHARIN
C
C     * AREA-WEIGHTED INTERPOLATION OF SCALAR OR VECTOR FIELDS
C     * FROM ONE GRID TO ANOTHER.
C
C     * UG,VG         = INPUT FIELDS.
C     * ILG1,ILG,ILAT = DIMENSIONS OF INPUT GRIDS.
C     *                 ILG1=ILG+1 IF THERE IS CYCLIC LONG.
C                       OTHERWISE ILG1=ILG.
C     * SLON          = LON (DEG) OF INPUT GRID COLUMNS.
C     * SLAT          = LAT (DEG) OF INPUT GRID ROWS.
C     * SLONH         = HALF LON (DEG) OF INPUT GRID COLUMNS
C                       (LON BOUNDARIES OF INPUT GRID CELLS)
C     * SLATH         = HALF LAT (DEG) OF INPUT GRID ROWS
C                       (LAT BOUNDARIES OF INPUT GRID CELLS)
C     * UL,VL         = OUTPUT FIELDS.
C     * NLG1,NLG,NLAT = DIMENSIONS OF OUTPUT GRID.
C     * DLON          = LON (DEG) OF OUTPUT GRID COLUMNS.
C     * DLAT          = LAT (DEG) OF OUTPUT GRID ROWS.
C     * DLONH         = HALF LON (DEG) OF OUTPUT GRID COLUMNS.
C                       (LON BOUNDARIES OF OUTPUT GRID CELLS)
C     * DLATH         = HALF LAT (DEG) OF OUTPUT GRID ROWS.
C                       (LAT BOUNDARIES OF OUTPUT GRID CELLS)
C     * WL            = AREA FRACTION WITH NON-MISSING VALUES.
C     * INTERP        = 5: AERA-WIGHTED INTERPOLATION.
C                       (THIS IS THE ONLY VALID OPTION.)
C     * IFLD          = 0 UG IS SCALAR FIELDS.
C                       (THIS IS THE ONLY VALID VALUE SO FAR)
C
C      IMPLICIT REAL (A-H,O-Z),
C     +INTEGER (I-N)
      IMPLICIT NONE
      REAL USP,VSP,UNP,VNP,U,V
      INTEGER ILG,ILG1,ILAT,NLG,NLG1,NLAT,INTERP,IFLD,I,J,
     1     I1,I2,J1,J2,II,JJ,II1,LUSP,LVSP,LUNP,LVNP
      REAL UG(ILG1,ILAT),VG(ILG1,ILAT),
     1     UL(NLG1,NLAT),VL(NLG1,NLAT),WL(NLG1,NLAT)
       REAL SLAT(ILAT),SLON(ILG1),DLAT(NLAT),DLON(NLG1),
     1     SLATH(ILAT+1),SLONH(ILG+1),DLATH(NLAT+1),DLONH(NLG+1)
      REAL R2D,D2R,SLN1,SLN2,DX,DY,DW,WTOT,WGHT,S

      REAL PI,SPVAL,SPVALT,SMALL
      DATA PI/3.14159265358979E0/,SPVAL/1.E+38/,
     1     SPVALT/1.E+32/,SMALL/1.E-12/
C-----------------------------------------------------------------------
      R2D=180.E0/PI
      D2R=PI/180.E0
      IF(INTERP.NE.5) CALL                         XIT('GGIGGW',-1)
      IF(IFLD.NE.0) THEN
        WRITE(6,'(2A)')
     1       '*** ERROR: AREA-WEIGHTED INTERPOLATION OF VECTOR FIELDS',
     2       ' IS NOT IMPLEMENTED YET. PLEASE USE BILINEAR OR CUBIC.'
        CALL                                       XIT('GGIGGW',-2)
      ENDIF

      DO J=1,NLAT
C
C       * FIND LATITUDES OF OVERLAPPING INPUT GRID BOXES
C
        DO J1=1,ILAT
          IF(SLATH(J1+1).GT.DLATH(J))GOTO 100
        ENDDO
 100    CONTINUE
        DO J2=ILAT,J1,-1
          IF(SLATH(J2).LT.DLATH(J+1))GOTO 110
        ENDDO
 110    CONTINUE

        DO I=1,NLG
C
C         * FIND LONGITUDES OF OVERLAPPING INPUT GRID BOXES
C
          DO I1=1-ILG,2*ILG
            IF(I1+1.LT.1)THEN
              SLN1=SLONH(I1+1+ILG)-360.E0
            ELSE IF(I1+1.GT.ILG)THEN
              SLN1=SLONH(I1+1-ILG)+360.E0
            ELSE
              SLN1=SLONH(I1+1)
            ENDIF
            IF(SLN1.GT.DLONH(I))GOTO 120
          ENDDO
 120      CONTINUE
          DO I2=2*ILG,I1,-1
            IF(I2.LT.1)THEN
              SLN2=SLONH(I2+ILG)-360.E0
            ELSE IF(I2.GT.ILG)THEN
              SLN2=SLONH(I2-ILG)+360.E0
            ELSE
              SLN2=SLONH(I2)
            ENDIF
            IF(SLN2.LT.DLONH(I+1))GOTO 130
          ENDDO
 130      CONTINUE
C
C         * CALCULATE AREA AVERAGE FOR POINT (I,J)
C
          WTOT=0.E0
          WGHT=0.E0
          S=0.E0
          DO JJ=J1,J2
            DY=SIND(MIN(SLATH(JJ+1),DLATH(J+1)))-
     1         SIND(MAX(SLATH(JJ  ),DLATH(J  )))
            DO II=I1,I2
              IF(II.LT.1)THEN
                SLN1=SLONH(II+ILG)-360.E0
                II1=II+ILG
              ELSE IF(II.GT.ILG)THEN
                SLN1=SLONH(II-ILG)+360.E0
                II1=II-ILG
              ELSE
                SLN1=SLONH(II)
                II1=II
              ENDIF
              IF(II+1.LT.1)THEN
                SLN2=SLONH(II+1+ILG)-360.E0
              ELSE IF(II.GT.ILG)THEN
                SLN2=SLONH(II+1-ILG)+360.E0
              ELSE
                SLN2=SLONH(II+1)
              ENDIF
              DX=MIN(SLN2,DLONH(I+1))-MAX(SLN1,DLONH(I))
              DW=DX*DY
              IF(ABS(UG(II1,JJ)-SPVAL).GT.SPVALT)THEN
                WGHT=WGHT+DW
                S=S+UG(II1,JJ)*DW
              ENDIF
              WTOT=WTOT+DW
            ENDDO
          ENDDO
          WL(I,J)=WGHT/WTOT
          IF(WL(I,J).LT.SMALL)THEN
            UL(I,J)=SPVAL
          ELSE
            UL(I,J)=S/WGHT
          ENDIF
        ENDDO
C
C       * CYCLIC LONGITUDE
C
        IF (NLG1.NE.NLG) THEN
          UL(NLG1,J)=UL(1,J)
          WL(NLG1,J)=WL(1,J)
        ENDIF
      ENDDO
      RETURN
      END
