SUBROUTINE FBARN(XLATN,NR,NC,XGA,XGD,XGE,XDA,
     *                  XROW,XCOL,NSTA,XINC,IFREQ,IGUESS)
C $ SUBROUTINE FBARN(XLATN,NR,NC,XGA,XGD,XGE,XDA,
C $ *XROW,XCOL,NSTA,XINC,IFREQ,IGUESS) (WLH)
C $ FAST VERSION OF BARNES INTERPOLATION OF RANDOMLY LOCATED OBSERVATIONS
C $   INTO UNIFORM LATITUDE-LONGITUDE GRID FOR ONE SET OF OBSERVATIONS.
C $ INPUT:
C $   XLATN = (R) NORTH LATITUDE OF GRID TOP IN DEGREES
C $   NR = (I) # ROWS IN GRID
C $   NC = (I) # COLUMNS IN GRID
C $   XDA =(R) ARRAY XDA(NSTA) CONTAINING OBSERVATIONS
C $   XROW = (R) ARRAY XROW(NSTA) CONTAINING ROW LOCATIONS OF OBSERVATIONS
C $   XCOL = (R) ARRAY XCOL(NSTA) CONTAINING COLUMN LOCATIONS OF OBSERVATIONS
C $      XROW & XCOL ARE MEASURED AS GRID ROW,COL . THUS AN OBSERVATION
C $      AT THE UPPER LEFT GRID POINT WOULD HAVE XROW & XCOL BOTH = 1.0
C $              ****** WARNING ******   THE VALUES IN XDA, XROW & XCOL
C $               ARE CHANGED BY THIS ROUTINE
C $   NSTA = (I) # OF OBSERVATIONS ( STATIONS )
C $   XINC = (R) LATITUDE BETWEEN ROWS, IN DEGREES
C $   IFREQ = (I) LOWPASS FILTER BANDWIDTH PARAMETER, INVERSELY
C $     PROPORTIONAL TO SPATIAL FREQUENCY.  GOOD VALUES FOR IFREQ DEPEND
C $     ON STATION DENSITY PER GRID SQUARE.  TOO LOW A VALUE MAY CAUSE A
C $     FLOATING POINT ABORT.  VALUES FROM 30 TO 60 ARE COMMON ON MCIDAS
C $   IGUESS = (I) NON-ZERO IF GUESS GRID INPUT THROUGH XGA
C $ OUTPUT:
C $   XGA = (R) ARRAY XGA(NR, NC), INTERPOLATED FROM XDA
C $ WORK:
C $   XGD = (R) ARRAY XGD(NR,NC) FOR WEIGHT SUM DENOMINATOR FOR 1ST PASS
C $   XGE = (R) ARRAY XGE(NR,NC) FOR WEIGHT SUM DENOMINATOR FOR 2ND PASS
C $$ FBARN = OBSERVATION, GRID, COMPUTATION
C
C   ALGORITHM AND SUBROUTINE BY BILL HIBBARD, SEPT 1979.
C THE COMPUTING TIME OF THIS ALGORITHM IS PROPORTIONAL TO
C NR*NC+NSTA*LOG(NSTA), COMPARED TO NR*NC*NSTA FOR INTPO.
C THE PARAMETERS ARE IDENTICAL TO THE PARAMETERS FOR THE INTPO ROUTINE
C IN THE MCIDAS SVCLIB LIBRARY, EXCEPT FOR THE ADDITION OF THE WORK
C ARRAYS XGD & XGE AND THE INPUT VARIABLE IGUESS.  IN PARTICULAR, IF
C IFREQ IS THE SAME IN THIS ROUTINE AND IN INTPO, THEN THEY SHOULD GIVE
C THE SAME INTERPOLATIONS, UP TO THE APPROXIMATION.
C
C-----*** WARNING *** ***(WARNING!)************
C THE VALUES IN XDA, XROW & XCOL ARE CHANGED BY THIS ROUTINE
C
      PARAMETER (MXROW=120,MXCOL=120)

      REAL*8 DEG2RAD, PI

      DIMENSION XGA(*),XGD(*),XGE(*)
      DIMENSION XDA(*)
      DIMENSION XROW(*),XCOL(*)
      COMMON/CONS/NRP
C THE DIMENSION OF IST MUST BE AT LEAST NR+2
      DIMENSION IST(MXROW+2)
      DIMENSION EXT(102)
C THE DIMENSION OF DBA,   DGA,   DBD, DBE, DGD, AND
C DGE MUST ALL BE AT LEAST NC
      DIMENSION DBA(MXCOL),DGA(MXCOL)
      DIMENSION DBD(MXCOL),DBE(MXCOL),DGD(MXCOL),DGE(MXCOL)
C THESE CONSTANTS DEFINE AN APPROXIMATION TO EXP(-X**2)
      DATA CLAMB/1.12045/,BETA/2.34/,GAMA/2.7495/
C THIS CONSTANT DEFINES THE CHANGE IN BANDWIDTH FROM THE FIRST TO SECOND
C PASS
      DATA SECND/1.825741858/
      PI = 4.D0*ATAN(1.D0)
      DEG2RAD = PI/180.D0
C THIS CONSTANT DEFINES THE INTEGER SCALING
      NRP=NR+1
      IF(NR .GT. MXROW .OR. NR .LT. 1) GO TO 999
      IF(NC .GT. MXCOL .OR. NC .LT. 1) GO TO 999
      IF(NSTA .LT. 2) GO TO 999
C
C
C TWO DIMENSIONAL HEAP SORT OF STATIONS. GROUPS BETWEEN HORIZONTAL
C GRID LINES ARE SORTED BY LONGITUDE. THE GROUPS ARE ORDERED BY
C LATITUDE
      NSTAH=NSTA/2
      DO 1 NN=1,NSTAH
      CALL HEAPFY(XDA,XROW,XCOL,NSTAH+1-NN,NSTA)
1     CONTINUE
      DO 2 NN=2,NSTA
      CALL INTER(XDA,XROW,XCOL,1,NSTA+2-NN)
      CALL HEAPFY(XDA,XROW,XCOL,1,NSTA+1-NN)
2     CONTINUE
C GENERATE POINTERS TO STARTS OF GROUPS OF STATIONS BETWEEN HORIZONTAL
C GRID LINES
      NSTAP=NSTA+1
      NRPP=NR+2
      DO 3 I=1,NRPP
3     IST(I)=NSTAP
      DO 4 N=1,NSTA
      K=XROW(N)+1.
      IF(K .LT. 1) K=1
      IF(K .GT. NRP) K=NRP
      IF(IST(K) .GT. N) IST(K)=N
4     CONTINUE
      MAX=NSTAP
      DO 5 II=1,NRP
      I=NRP+1-II
      IF(IST(I) .LE. MAX) GO TO 5
      IST(I)=MAX
5     MAX=IST(I)
C CALCULATE EXPONENTS FROM LOW PASS BANDWIDTH PARAMETER IFREQ
      FREQ=IFREQ
      FF=SQRT(10.0/FREQ)
      FFA=FF*SECND
      SECI=1.0/SECND
      BG=BETA/GAMA
      BC=FF*BETA
      GC=FF*GAMA
      BCA=FFA*BETA
      GCA=FFA*GAMA
      BCF=SECI*BG
      GCF=SECI
      BCAF=BG
      GCAF=1.0
C FILL LOOKUP TABLE FOR SMALL ARGUMENT EXPONENTIALS
      ECN=-GCA/100.0
      DO 7 I=1,102
      VAL=(I-1)*ECN
      EXT(I)=EXP(VAL)
7     CONTINUE
C ZERO ARRAYS FOR GRIDDED DATA
      NRNC=NR*NC
      DO 6 I=1,NRNC
C***GUESS MOD
      IF(IGUESS.EQ.0) XGA(I)=0.0
C***
      XGD(I)=0.0
6     XGE(I)=0.0
C GET MULTIPLY COEFFICIENTS FROM EXPONENTS. THESE ARE FOR FILTER ALONG
C LINES OF LONGITUDE ( COLUMNS )
      EBC=EXP(-BC)
      EGC=EXP(-GC)
      EBCA=EXP(-BCA)
      EGCA=EXP(-GCA)
C TWO PASSES, 1-INITIAL INTERPOLATION, 2-CORRECTION
      DO 100 IPASS=1,2
C SCAN DATA, 1-DECREASING LATITUDE, 2- INCREASING LATITUDE
      DO 99 IUD=1,2
      IF(IUD .EQ. 2) GO TO 8
      JLO=0
      JINC=1
      FLAT=(XLATN+XINC)
      FLINC=-XINC
      GO TO 9
8     JLO=NR+1
      JINC=-1
      FLAT=(XLATN-XINC*NR)
      FLINC=XINC
C *****                                          *****
9     CONTINUE
C ZERO RUNNING VERTICAL SUMS
      DO 10 I=1,NC
      DBA(I)=0.0
      DGA(I)=0.0
10    CONTINUE
      IF(IPASS .EQ. 2) GO TO 12
      DO 11 I=1,NC
      DBD(I)=0.0
      DBE(I)=0.0
      DGD(I)=0.0
11    DGE(I)=0.0
12    CONTINUE
      J=JLO
C VERTICAL SCAN LOOP
      DO 98 JJ=1,NR
      J=J+JINC
      JHUN=J*100
C ADJUST EXPONENTS AND MULTIPLIERS FOR LATITUDE DEPENDENT DISTANCE
C DISTORTIONS. THESE ARE FOR FILTERING ALONG LINES OF LATITUDE ( ROWS )
      FLAT=FLAT+FLINC
C
C--- PREVENT ROUND OFF ERROR AND ASSURE THAT COS(90) = 0
C
      IF(ABS(FLAT) .EQ. 90.) THEN
        CL = 0.0
      ELSE
        CL=COS(FLAT*DEG2RAD)
      ENDIF
      IF(CL .LT. 0.0) GO TO 999
      BR=BC*CL
      GR=GC*CL
      BRA=BCA*CL
      GRA=GCA*CL
      BRF=BCF*CL
      GRF=GCF*CL
      BRAF=BCAF*CL
      GRAF=GCAF*CL
      EBR=EXP(-BR)
      EGR=EXP(-GR)
      EBRA=EXP(-BRA)
      EGRA=EXP(-GRA)
C DECAY VERTICAL RUNNING SUMS
      DO 13 I=1,NC
      DBA(I)=EBC*DBA(I)
      DGA(I)=EGC*DGA(I)
13    CONTINUE
      IF(IPASS .EQ. 2) GO TO 15
      DO 14 I=1,NC
      DBD(I)=EBC*DBD(I)
      DGD(I)=EGC*DGD(I)
      DBE(I)=EBCA*DBE(I)
14    DGE(I)=EGCA*DGE(I)
15    CONTINUE
C SCAN DATA, 1-INCREASING LONGITUDE, 2-DECREASING LONGITUDE
      DO 97 IRL=1,2
      IF(IRL .EQ. 2) GO TO 16
      ILO=0
      IINC=1
      GO TO 17
16    ILO=NC+1
      IINC=-1
17    CONTINUE
C ZERO RUNNING HORIZONTAL SUMS
      HBBA=0.0
      HBGA=0.0
      HGBA=0.0
      HGGA=0.0
      IF(IPASS .EQ. 2) GO TO 18
      HBBD=0.0
      HBBE=0.0
      HBGD=0.0
      HBGE=0.0
      HGBD=0.0
      HGBE=0.0
      HGGD=0.0
      HGGE=0.0
18    CONTINUE
C INITIALIZE POINTERS INTO STATION GROUP AT THIS LATITUDE
      IS=J+IUD-1
      IR=3-2*IRL
      IC=3-2*IUD
      KSTAR=J+IUD+IRL-2
      LSTAR=J+IUD-IRL+1
      KSTA=IST(KSTAR)+1-IRL
      LSTA=IST(LSTAR)-2+IRL
      IJ=J+NR*(ILO-1)
      IJINC=IINC*NR
      I=ILO
C HORIZONTAL SCAN LOOP
      DO 96 II=1,NC
      I=I+IINC
      IHUN=I*100
      IJ=IJ+IJINC
C DECAY HORIZONTAL SUMS
      HBBA=EBR*HBBA
      HBGA=EGR*HBGA
      HGBA=EBR*HGBA
      HGGA=EGR*HGGA
      IF(IPASS .EQ. 2) GO TO 19
      HBBD=EBR*HBBD
      HBBE=EBRA*HBBE
      HBGD=EGR*HBGD
      HBGE=EGRA*HBGE
      HGBD=EBR*HGBD
      HGBE=EBRA*HGBE
      HGGD=EGR*HGGD
      HGGE=EGRA*HGGE
19    CONTINUE
C START OF LOOP TO FIND ALL STATIONS IN THIS GRID BOX, ADJUST
C OBSERVATION VALUES BY EXPONENT LOOK UP TABLE, AND ADD VALUES TO
C HORIZONTAL RUNNING SUMS
24    IF(IR*(KSTA-LSTA) .GT. 0) GO TO 25
      NU=IR*(IHUN-100*XCOL(KSTA))
      IF(NU .LT. 0) GO TO 25
      MV=IC*(JHUN-100*XROW(KSTA))
      IF(MV .GT. 0) GO TO 40
      KSTA=KSTA+IR
      GO TO 24
40    CONTINUE
C SCALE FOR DISTANCE FROM STATION TO GRID POINT
      CBCMV=BCF*MV
      CGCMV=GCF*MV
      CBRNU=BRF*NU
      CGRNU=GRF*NU
C CHECK IF LOOKUP TABLE CAN BE USED
      IF(MV .GE. 100) GO TO 20
      IF(NU .GE. 100) GO TO 20
      CGBC=CGCMV-CBCMV
      CGBR=CGRNU-CBRNU
      IBC=CBCMV+1.5
      IGCC=CGBC+1.5
      IBR=CBRNU+1.5
      IGR=CGBR+1.5
      RBC=EXT(IBC)
      RGC=EXT(IGCC)*RBC
      RBR=EXT(IBR)
      RGR=EXT(IGR)*RBR
      RBRCL=RBR*CLAMB
      RBB=RBC*RBRCL
      RBG=RBC*RGR
      RGB=RGC*RBRCL
      RGG=RGC*RGR
      GO TO 21
20    RBB=(CBCMV+CBRNU)*ECN
      RBG=(CBCMV+CGRNU)*ECN
      RGB=(CGCMV+CBRNU)*ECN
      RGG=(CGCMV+CGRNU)*ECN
      RBB=EXP(RBB)*CLAMB
      RBG=EXP(RBG)
      RGB=EXP(RGB)*CLAMB
      RGG=EXP(RGG)
21    SA=XDA(KSTA)
      KSTA=KSTA+IR
C ADD VALUES TO HORIZONTAL RUNNING SUMS
      HBBA=HBBA+RBB*SA
      HBGA=HBGA+RBG*SA
      HGBA=HGBA+RGB*SA
      HGGA=HGGA+RGG*SA
      IF(IPASS .EQ. 2) GO TO 24
      HBBD=HBBD+RBB
      HBGD=HBGD+RBG
      HGBD=HGBD+RGB
      HGGD=HGGD+RGG
      CBCMV=BCAF*MV
      CGCMV=GCAF*MV
      CBRNU=BRAF*NU
      CGRNU=GRAF*NU
      IF(MV .GE. 100) GO TO 22
      IF(NU .GE. 100) GO TO 22
      CGBC=CGCMV-CBCMV
      CGBR=CGRNU-CBRNU
      IBC=CBCMV+1.5
      IGCC=CGBC+1.5
      IBR=CBRNU+1.5
      IGR=CGBR+1.5
      RBC=EXT(IBC)
      RGC=EXT(IGCC)*RBC
      RBR=EXT(IBR)
      RGR=EXT(IGR)*RBR
      RBRCL=RBR*CLAMB
      RBB=RBC*RBRCL
      RBG=RBC*RGR
      RGB=RGC*RBRCL
      RGG=RGC*RGR
      GO TO 23
22    RBB=(CBCMV+CBRNU)*ECN
      RBG=(CBCMV+CGRNU)*ECN
      RGB=(CGCMV+CBRNU)*ECN
      RGG=(CGCMV+CGRNU)*ECN
      RBB=EXP(RBB)*CLAMB
      RBG=EXP(RBG)
      RGB=EXP(RGB)*CLAMB
      RGG=EXP(RGG)
23    HBBE=HBBE+RBB
      HBGE=HBGE+RBG
      HGBE=HGBE+RGB
      HGGE=HGGE+RGG
      GO TO 24
C END OF LOOP FOR STATIONS IN A GRID BOX
25    CONTINUE
C ADD HORIZONTAL SUMS TO RUNNING VERTICAL SUMS
      DBA(I)=DBA(I)+CLAMB*(HBBA-HBGA)
      DGA(I)=DGA(I)+HGBA-HGGA
      IF(IPASS .EQ. 2) GO TO 26
      DBD(I)=DBD(I)+CLAMB*(HBBD-HBGD)
      DBE(I)=DBE(I)+CLAMB*(HBBE-HBGE)
      DGD(I)=DGD(I)+HGBD-HGGD
      DGE(I)=DGE(I)+HGBE-HGGE
26    CONTINUE
96    CONTINUE
97    CONTINUE
C ADD VERTICAL SUMS TO GRID ARRAYS
      IJ=J-NR
      IF(IPASS .EQ. 2) GO TO 60
      DO 61 I=1,NC
      IJ=IJ+NR
      XGD(IJ)=XGD(IJ)+(DBD(I)-DGD(I))
      XGE(IJ)=XGE(IJ)+(DBE(I)-DGE(I))
61    CONTINUE
C***GUESS MOD
      IF(IGUESS.NE.0) GO TO 63
      IJ=J-NR
      DO 71 I=1,NC
      IJ=IJ+NR
      XGA(IJ)=XGA(IJ)+DBA(I)-DGA(I)
71    CONTINUE
C***
      GO TO 63
60    CONTINUE
      DO 62 I=1,NC
      IJ=IJ+NR
      IF(XGE(IJ).GT.0.0) THEN
         EFACT=1./XGE(IJ)
         XGA(IJ)=XGA(IJ)+(DBA(I)-DGA(I))*EFACT
      ELSE
         XGA(IJ)=1.E35
      ENDIF
62    CONTINUE
63    CONTINUE
98    CONTINUE
99    CONTINUE
      IF(IPASS .EQ. 2) GO TO 100
C ADJUST STATION OBSERVATIONS AFTER FIRST PASS FOR ERROR ADJUSTMENT
C CALCULATION ON SECOND PASS
      BC=BCA
      GC=GCA
      BCF=BCAF
      GCF=GCAF
      EBC=EBCA
      EGC=EGCA
      DO 28 I=1,NRNC
      IF(XGD(I).LE.1.E-20) XGD(I)=-99.
      IF(XGE(I).LE.1.E-20) XGE(I)=-99.
28    CONTINUE
C***GUESS MOD
      IF(IGUESS.NE.0) GO TO 250
      DO 29 IJ=1,NRNC
      IF(XGD(IJ).GT.0) THEN
         EFACT=1./XGD(IJ)
         XGA(IJ)=XGA(IJ)*EFACT
      ELSE
         XGA(IJ)=1.E35
      ENDIF
29    CONTINUE
250   CONTINUE
      IL=100
      IH=100*NC
      JL=100
      JH=100*NR
      DO 30 N=1,NSTA
      I=NINT(XCOL(N))
      J=NINT(XROW(N))
      IF(I .LT. 1) I=1
      IF(I .GT. NC) I=NC
      IF(J .LT. 1) J=1
      IF(J .GT. NR) J=NR
      IJ=J+NR*(I-1)
      XDA(N)=XDA(N)-XGA(IJ)
30    CONTINUE
100   CONTINUE
      RETURN
999   CONTINUE
C **** output error message here for: 'FBARN CALLED WITH BAD PARAMETERS' ****
      CALL ABORT(0)
      RETURN
      END
C
      SUBROUTINE HEAPFY(XDA,XROW,XCOL,LL,K)
C $ CALL HEAPFY(XDA, XROW, XCOL, LL, K)  (WLH)
C $ REORDER ELEMENTS OF HEAP AS PART OF HEAP SORT
C $ XDA =(R) ARRAY IDA(NSTA) CONTAINING FIRST OBSERVED PARAMETER
C $ XROW = (R) ARRAY XROW(NCOL) CONTAINING LATITUDES OF OBSERVATIONS
C $ XCOL = (R) ARRAY XCOL(NSTA) CONTAINING LONGITUDES OF OBSERVATIONS
C $ LL = (I) INPUT  STARTING LOCATION IN HEAP FOR REORDERING
C $ K = (I) INPUT  ENDING LOCATION IN HEAP FOR REORDER
C $$ HEAPFY = OBSERVATION, GRID, COMPUTATION
      DIMENSION XDA(*),XROW(*),XCOL(*)
      IFLAG=0
      I=LL
4     M=I + I
      IF(M .GT. K) RETURN
      N=M+1
      IF(ICOMP(XROW,XCOL,I,M) .EQ. 1) GO TO 2
      IF(N .GT. K) GO TO 3
      IF(ICOMP(XROW,XCOL,M,N) .EQ. 1) GO TO 3
      J=N
      GO TO 1
3     J=M
      GO TO 1
2     IF(N .GT. K) RETURN
      IF(ICOMP(XROW,XCOL,I,N) .EQ. 1) RETURN
      J=N
      GO TO 1
      ENTRY INTER(XDA,XROW,XCOL,LL,K)
C $ ENTRY INTER(XDA, XROW, XCOL, LL, K)  (WLH)
C $ INTERCHANGES (SWAPS) POSITIONS OF VECTORS
C $ XDA =(R) ARRAY IDA(NSTA) CONTAINING FIRST OBSERVED PARAMETER
C $ XROW = (R) ARRAY XROW(NCOL) CONTAINING LATITUDES OF OBSERVATIONS
C $ XCOL = (R) ARRAY XCOL(NSTA) CONTAINING LONGITUDES OF OBSERVATIONS
C $ LL = (I) INPUT  STARTING LOCATION IN HEAP FOR REORDERING
C $ K = (I) INPUT  ENDING LOCATION IN HEAP FOR REORDER
C $$ INTER = OBSERVATION, GRID, COMPUTATION
      IFLAG=1
      I=LL
      J=K
1     X=XDA(I)
      XDA(I)=XDA(J)
      XDA(J)=X
      X=XROW(I)
      XROW(I)=XROW(J)
      XROW(J)=X
      X=XCOL(I)
      XCOL(I)=XCOL(J)
      XCOL(J)=X
      IF(IFLAG .EQ. 1) RETURN
      I=J
      GO TO 4
      END
C
      FUNCTION ICOMP(XROW,XCOL,I,J)
C $ ICOMP(XROW, XCOL, I, J)  (WLH)
C $ PERFORMS LEXICAL COMPARISON OF 2 (LAT,LON) VECTORS
C $ XROW = (R) ARRAY XROW(NCOL) CONTAINING LATITUDES OF OBSERVATIONS
C $ XCOL = (R) ARRAY XCOL(NSTA) CONTAINING LONGITUDES OF OBSERVATIONS
C $ I = (I) INPUT INDEX OF FIRST VECTOR FOR COMPARE
C $ J = (I) INPUT INDEX OF SECOND VECTOR FOR COMPARE
C $$ ICOMP = OBSERVATION, GRID, COMPUTATION
      COMMON/CONS/NRP
      DIMENSION XROW(*),XCOL(*)
      K=XROW(I)+1
      IF(K .LT. 1) K=1
      IF(K .GT. NRP) K=NRP
      L=XROW(J)+1
      IF(L .LT. 1) L=1
      IF(L .GT. NRP) L=NRP
      IF(K .LT. L) GO TO 1
      IF(K .GT. L) GO TO 2
      IF(XCOL(I) .LE. XCOL(J)) GO TO 1
2     ICOMP=1
      RETURN
1     ICOMP=0
      RETURN
      END