SUBROUTINE CLOS2(KFILDO,CCALL,XP,YP,LTAG,ELEV,XLAPSE,
1 LNDSEA,XDATA,NSTA,K,R,LIO,
2 LS1,LS2,DIST1,DIST2,ELD,
3 ELDS,ELDXCL,IER)
C
C NOVEMBER 2007 GLAHN MDL
C MODIFIED FROM LAMP CLOS MAY 1993
C DECEMBER 2007 GLAHN MODIFIED FOR LIO AND DIFFERENT
C DATA TYPES
C DECEMBER 2007 GLAHN DEFINED AND USED RSQ; INITIALIZED
C LS1 AND LS2
C DECEMBER 2007 GLAHN MODIFIED TESTS FOR LIO ABOVE
C DISTSQ= (12/26/07)
C DECEMBER 2007 GLAHN ADDED ELD, ELDS, ELDXCL TO OUTPUT,
C ELEV TO INPUT.
C DECEMBER 2007 GLAHN ADDED 190 AND GO TO 190
C JANUARY 2008 GLAHN CHANGED USE OF ELDS IN CALL; ADDED
C XLAPSE
C FEBRUARY 2008 GLAHN CHANGED INDICES IN ELDS= BELOW 170
C MARCH 2008 COSGROVE ADDED COMMA TO FORMAT 170 FOR IBM
C FEBRUARY 2009 GLAHN CHANGED LNDSEA(K) TO LNDSEA(L) CHECK
C DECEMBER 2009 GLAHN IMPROVED COMMENT FOR ELDS
C
C PURPOSE
C TO FIND 2 CLOSEST NEIGHBORS OF STATION OF LOCATION K
C IN CCALL( ) LIST FOR U155 AND TO RETURN THE POSITIONS
C IN THE LIST AND THE DISTANCES. LIO CAN BE OF TYPE 0 OR 9.
C FOR TYPE LIO = 0, USE STATION TYPES 0 AND 3;
C = 9, USE STATION TYPES 6 AND 9.
C ALSO COMPUTES ELD, ELDS, AND ELDXCL (SEE DEFINITIONS BELOW)
C
C DATA SET USE
C KFILDO - UNIT NUMBER FOR OUTPUT (PRINT) FILE. (OUTPUT)
C
C VARIABLES
C INPUT
C KFILDO = UNIT NUMBER FOR OUTPUT (PRINT) FILE. (INPUT)
C CCALL(K) = CALL LETTERS OF STATIONS (K=1,NSTA).
C (CHARACTER*8) (INPUT)
C XP(J) = X-POSITION OF STATION J ON GRID (J=1,NSTA).
C YP(J) = Y-POSITION OF STATION J ON GRID (J=1,NSTA).
C LTAG(J) = DON'T USE THIS STATION IF LTAG( ) GT 0
C (J=1,NSTA).
C ELEV(K) = ELEVATION OF STATIONS (K=1,NSTA). (INPUT)
C XLAPSE(K) = CALCULATED LAPSE RATE IN UNITS OF THE VARIABLE
C BEING ANALYZED PER M. (K=1,KSTA). (INPUT)
C LNDSEA(K) = LAND/SEA INFLUENCE FLAG FOR EACH STATION
C (K=1,ND1).
C 0 = WILL BE USED FOR ONLY OCEAN WATER (=0)
C GRIDPOINTS.
C 3 = WILL BE USED FOR ONLY INLAND WATER (=3)
C GRIDPOINTS.
C 6 = WILL BE USED FOR BOTH INLAND WATER (=3)
C AND LAND (=9) GRIDPOINTS.
C 9 = WILL BE USED FOR ONLY LAND (=9) GRIDPOINTS.
C (INPUT)
C XDATA(J) = THE DATA FOR WHICH TO FIND THE DENSITY
C (J=1,ND1).
C NSTA = NUMBER OF STATIONS IN LIST.
C K = POSITION OF STATION IN LIST FOR WHICH CLOSEST
C STATIONS ARE NEEDED.
C R = RADIUS OVER WHICH TO SEARCH. THIS IS RCLOS2
C IN CALLING PROGRAM. IF IT TURNS OUT, WATER
C STATIONS CANNOT BE FOUND IN THIS RADIUS,
C RCLOS2 CAN BE INCREASED IN U405A. THIS WOULD
C MEAN A LARGER SEARCH RADIUS FOR ALL TYPES
C OF STATIONS. IT WOULD ALSO MEAN IDS OF SOME
C VARIABLES WOULD CHANGE.
C LIO = THE VARIABLE TO SPECIFY WHETHER THIS RUN IS
C FOR LAND (=9) OR WATER (=0) POINTS.
C
C OUTPUT
C LS1 = POSITION IN LIST OF CLOSEST STATION.
C LS2 = POSITION IN LIST OF 2ND CLOSEST STATION.
C DIST1 = DISTANCE FROM STATION K TO CLOSEST STATION.
C DIST2 = DISTANCE FROM STATION K TO 2ND CLOSEST STATION.
C ELD = ABSOLUTE DIFFERENCE IN ELEVATION BETWEEN
C THE STATION K AND THE CLOSEST STATION LS1
C ELDS = ABSOLUTE DIFFERENCE BETWEEN STATION DATA
C VALUE AND VALUE OF NEIGHBOR AFTER LAPSE
C AT THE CLOSEST STATION IS APPLIED.
C ELDXCL = THE PRODUCT OF THE DISTANCE TO THE CLOSEST
C STATION AND THE ELEVATION DIFFERENCE.
C IER = STATUS RETURN
C 0 = GOOD.
C 777 = CAN'T FIND TWO CLOSEST STATIONS WITHIN
C DISTANCE R.
C
C INTERNAL
C SAV1 = DISTANCE SQUARED TO CLOSEST STATION.
C SAV2 = DISTANCE SQUARED TO 2ND CLOSEST STATION.
C DISTSQ = DISTANCE (IN GRID UNITS) SQUARED BETWEEN
C TWO STATIONS.
C 1 2 3 4 5 6 7 X
C
C NONSYSTEM SUBROUTINES CALLED
C NONE.
C
CHARACTER*8 CCALL(NSTA)
C
DIMENSION XP(NSTA),YP(NSTA),LTAG(NSTA),LNDSEA(NSTA),XDATA(NSTA),
1 ELEV(NSTA),XLAPSE(NSTA)
C
CD CALL TIMPR(KFILDO,KFILDO,'START CLOS2 ')
C
IER=0
SAV1=999999.
SAV2=999999.
LS1=999999
LS2=999999
RSQ=R*R+.01
C THE SMALL CONSTANT IS ADDED TO ASSURE A POINT IS NOT
C ELIMINATED BECAUSE OF ROUNDOFF.
C
DO 162 L=1,NSTA
IF(XDATA(L).GT.9998.5)GO TO 162
C DO NOT USE STATION WITH MISSING DATA.
C
IF(K.EQ.L)GO TO 162
C DO NOT USE STATION ITSELF.
C
IF(ABS(XP(K)-XP(L)).GT.R)GO TO 162
IF(ABS(YP(K)-YP(L)).GT.R)GO TO 162
C IN A LONG LIST OF STATIONS, THE ABOVE TWO TESTS SHOULD BE
C MORE EFFICIENT THAT ALWAYS CALCULATING THE DISTANCE.
C ALSO, THEY SHOULD RULE OUT MORE THAN THE TWO FOLLOWING
C TESTS.
C
IF(LTAG(L).GT.0)GO TO 162
C
IF(LIO.EQ.9.AND.LNDSEA(L).LT.6)GO TO 162
IF(LIO.EQ.0.AND.LNDSEA(L).GT.3)GO TO 162
C THESE TESTS USE LAND STATIONS OR JOINT LAND/INLAND WATER
C STATIONS FOR LIO=9 (TO BE USED FOR LAND GRIDPOINTS),
C AND BOTH OCEAN AND INLAND WATER WHEN LIO = 0 (TO BE USED
C FOR ALL WATER GRIDPOINTS). UPON ENTRY, THE STATION K
C IS OF THE CORRECT TYPE; THIS CHECK IS ON STATION L.
C
DISTSQ=(XP(K)-XP(L))**2+(YP(K)-YP(L))**2
C
IF(DISTSQ.LT.RSQ)THEN
IF(SAV2.LE.DISTSQ)GO TO 162
IF(SAV1.GT.DISTSQ)GO TO 1611
SAV2=DISTSQ
LS2=L
GO TO 162
C
1611 SAV2=SAV1
SAV1=DISTSQ
LS2=LS1
LS1=L
ENDIF
C
162 CONTINUE
C
IF(SAV2.GT.999998.5)THEN
IER=777
WRITE(KFILDO,170)CCALL(K)
170 FORMAT(/,' ****CANNOT FIND TWO CLOSEST STATIONS TO =',A8,
1 ' IN CLOS2.')
GO TO 190
ELSE
DIST1=SQRT(SAV1)
DIST2=SQRT(SAV2)
C
CD WRITE(KFILDO,1615)LS1,LS2,CCALL(LS1),CCALL(LS2),DIST1,DIST2
CD1615 FORMAT(/' AT 1615 IN CLOS2--',
CD 1 'LS1,LS2,CCALL(LS1),CCALL(LS2),SAV1,SAV2',2I6,2A8,2F10.2)
ENDIF
C
ELD=ABS(ELEV(K)-ELEV(LS1))
ELDXCL=ELD*DIST1
ELDS=ABS(XDATA(K)-(XDATA(LS1)+XLAPSE(LS1)*(ELEV(K)-ELEV(LS1))))
C ELDS IS THE ABSOLUTE DIFFERENCE BETWEEN THE DATA VALUE AT K
C AND THE VALUE ESTIMATED FROM THE CLOSEST STATION, CONSIDERING
C THE LAPSE RATE AND THE DIFFERENCES IN ELEVATION. (IN
C EVALUATION IN CLOS2G, THE LOCATION K WILL BE THE GRIDPOINT
C AND THE CLOSEST STATION WILL BE USED TO FIND THE DIFFERENCE
C BETWEEN THE ANALYSIS VALUE AND THE ESTIMATED VALUE.)
C
CD WRITE(KFILDO,180)CCALL(K),CCALL(LS1),CCALL(LS2),DIST1,DIST2
CD180 FORMAT(/,' STATION AND 2 CLOSEST ARE',3(3X,A8),
CD 1 ' AND DISTANCES ARE',2F6.1,' GRID UNITS.')
CD WRITE(KFILDO,181)ELD,ELDXCL,XDATA(K),XDATA(LS1),
CD 1 ELEV(K),ELEV(LS1),XLAPSE(K)
CD181 FORMAT(' IN CLOS2 AT 181--ELD,ELDXCL,XDATA(K),XDATA(LS1),',
CD 1 'ELEV(K),ELEV(LS1),XLAPSE(K)',7F7.1)
C
CD CALL TIMPR(KFILDO,KFILDO,'END CLOS2 ')
190 RETURN
END