SUBROUTINE VARILG(KFILDO,ID,
1 CCALL,XP,YP,LTAG,LNDSEA,XDATA,ELEV,XLAPSE,NSTA,
2 ERROR,NX,NY,
3 CPNDFD,SEALND,TELEV,NXE,NYE,
2 ITABLE,RTABLE,IDIM,JDIM,M,RELVAR,VARTAB,
5 ITERML,ITERMW,RMULT,ISTOP,IER)
C
C FEBRUARY 2008 GLAHN MDL
C MARCH 2008 COSGROVE ADDED COMMAS TO 163,173 FOR IBM
C FEBRUARY 2009 GLAHN ADDED TO PURPOSE; CHANGED
C IF(LNDSEA(K).LT.6) TO GT AND
C INCREASED RSQW + .01 AT 170
C MARCH 2009 GLAHN ADDED VARI=0 ABOVE DO 182
C MARCH 2009 GLAHN CHANGED IF(LTAG(K).GT.0) TO
C IF(LTAG(K).NE.0) IN TWO PLACES
C MARCH 2009 GLAHN CHANGED 2ND WORD STRUCTURE AND
C OTHER STUFF
C APRIL 2009 GLAHN INSERT MISSING WHEN 2 POINTS NOT
C FOUND AND PROCEED; ADDED ICOUNT
C JULY 2010 SCALLION OPEN MP STATEMENTS PUT IN IN
C 2 PLACES; SPELL CHECK; COMMENTS;
C ELIMINATED "SAVE IFIRST"
C SEPTEMBER 2011 IM/GLAHN INSERTED GO TO 199 IN DO 1635
C AND 1735 LOOPS
C NOVEMBER 2011 IM/GLAHN ADDED RMULT TO CALL;
C CORRECTED COMMENTS; ADDED IB=2
C PROCEDURES
C
C PURPOSE
C TO ACCUMULATE THE ANALYSIS ERROR IN ERROR( , ) FOR
C THE TERMS IN THE ERROR ESTIMATION EQUATIONS INVOLVING
C DATA DIFFERENCES IT DEALS SEPARATELY WITH BOTH LAND
C AND WATER EQUATIONS AND GRIDPOINTS. THE ERROR( , )
C IS ACCUMULATED FOR ALL TERMS INVOLVING 1397XXXX, 1497XXXX,
C 1597XXXX, OR 1697XXXX, WHERE XXXX IS A RADIUS R:
C
C 13970055 = AVERAGE OF THE ABSOLUTE VALUES OF THE
C DIFFERENCES BETWEEN THE DATA VALUES AND THEIR
C MEAN WITHIN 55 GRIDLENGTHS OF THE GRIDPOINT
C AFTER APPLYING THE LAPSE RATE, THE DIFFERENCE
C IN ELEVATION, AND QUADRATIC DISTANCE
C WEIGHTING,
C 14970045 = AVERAGE OF THE ABSOLUTE VALUES OF THE
C DIFFERENCES BETWEEN THE DATA VALUES AND THEIR
C MEAN WITHIN 45 GRIDLENGTHS OF THE GRIDPOINT
C AFTER APPLYING THE LAPSE RATE, THE DIFFERENCE
C IN ELEVATION, AND QUADRATIC DISTANCE WEIGHTING,
C 15970035 = AVERAGE OF THE ABSOLUTE VALUES OF THE
C DIFFERENCES BETWEEN THE DATA VALUES AND THEIR
C MEAN WITHIN 35 GRIDLENGTHS OF THE GRIDPOINT
C AFTER APPLYING THE LAPSE RATE, THE DIFFERENCE
C IN ELEVATION, AND QUADRATIC DISTANCE
C WEIGHTING, ANDNNNNNNNN
C 16970027 = AVERAGE OF THE ABSOLUTE VALUES OF THE
C DIFFERENCES BETWEEN THE DATA VALUES AND THEIR
C MEAN WITHIN 27 GRIDLENGTHS OF THE GRIDPOINT
C AFTER APPLYING THE LAPSE RATE, THE DIFFERENCE
C IN ELEVATION, AND QUADRADIC DISTANCE
C WEIGHTING.
C
C THE VALUES 55, 45, 35, AND 27 ARE EXAMPLES AND COULD BE
C DIFFERENT.
C
C THE CALCULATION IS MADE IN THE FOLLOWING WAY:
C 1) EACH STATION WITHIN THE RADIUS R OF THE GRIDPOINT
C IS ADJUSTED TO THE GRIDPOINT WITH ITS LAPSE RATE
C AND ITS VERTICAL DIFFERENCE FROM THE GRIDPOINT.
C 2) FOR EACH STATION, THE QUADRATIC DISTANCE WEIGHTING
C FUNCTION IS CALCULATED AND STORED IN TEMPW( ).
C 3) THE VALUES IN 1) ARE AVERAGED.
C 4) THE AVERAGE DIFFERENCE OF THE INDIVIDUAL VALUES
C COMPUTED IN 1) WEIGHTED BY THE WEIGHTING FUNCTION
C COMPUTED IN 2) AND THE AVERAGE COMPUTED IN 3)
C IS COMPUTED.
C
C THIS IS THE VARIABILITY OF THE DATA VALUES AFTER
C BEING ADJUSTED TO THE ELEVATION OF THE GRIDPOINT AND
C THE DISTANCE BY THE WEIGHTING FUNCTION.
C
C NOTE THAT WHEN THERE IS ONLY ONE STATION OR NONE WITHIN
C THE RADIUS, A BACKUP RADIUS IS USED FROM
C RELVAR(JJ)*VARTAB(L). ONLY WHEN THE LARGEST VALUE
C CANNOT YIELD TWO CLOSEST STATIONS DOES A POINT GET SET
C TO MISSING. THIS DOES NOT CAUSE AN ERROR RETURN.
C
C THE RADIUS OF SEARCH IS TAKEN FROM THE 2ND WORD OF THE
C EQUATION OR A BACKUP IS USED. THE LARGEST BACKUP VALUE
C SHOULD USUALLY BE LARGE ENOUGH UNLESS THERE ARE
C SIGNIFICANT MISSING DATA.
C
C OPEN MP ISTOP AND IER ARE NOT IN THE OPEN MP LOOP.
C IFIRST COULD BE CONFLICTED, BUT WOULD NOT HURT OPERATION.
C IT COULD CAUSE TWO OR MORE DIAGNOSTICS TO PRINT
C INSTEAD OF ONE.
C ICOUNT AND JCOUNT COULD BE CONFLICTED, BUT THEY ONLY
C CONTROL DIAGNOSTICS AT THE END AND IT IS NOT
C CRUCIAL THEY BE EXACT.
C
C DATA SET USE
C KFILDO - UNIT NUMBER FOR OUTPUT (PRINT) FILE. (OUTPUT)
C
C VARIABLES
C KFILDO = UNIT NUMBER FOR OUTPUT (PRINT) FILE. (INPUT)
C ID(J) = THE VARIABLE ID'S OF THE ANALYSIS (J=1,4).
C (INPUT)
C CCALL(K) = CALL LETTERS OF STATIONS (K=1,NSTA). USED
C ONLY FOR DIAGNOSTICS. (CHARACTER*8) (INPUT)
C XP(K) = X-POSITION OF STATION K ON GRID (J=1,NSTA).
C (INPUT)
C YP(K) = Y-POSITION OF STATION K ON GRID (J=1,NSTA).
C (INPUT)
C LTAG(K) = DON'T USE THIS STATION IF LTAG( ) NE 0
C (K=1,NSTA). STATIONS WITH DATA BUT TOSSED
C HAVE LTAG( ) = -1. POSSIBILITY OF LTAG( )
C = -3 NOT IMPLEMENTED. (INPUT)
C LNDSEA(K) = LAND/SEA INFLUENCE FLAG FOR EACH STATION
C (K=1,NSTA).
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. THESE ARE
C VERY RARE AND ARE LUMPED WITH WATER.
C 9 = WILL BE USED FOR ONLY LAND (=9) GRIDPOINTS.
C NOTE: ALL WATER STATIONS, INCLUDING
C LNDSEA( ) = 6, ARE GROUPED INTO ONE WATER
C CATEGORY.
C (INPUT)
C XDATA(K) = THE DATA FOR WHICH TO FIND THE VARIABILITY
C (K=1,NSTA). (INPUT)
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 NSTA = NUMBER OF STATIONS IN LIST. (INPUT)
C ERROR(IX,JY) = HOLDS THE ACCUMULATED ERROR (IX=1,NY) (JY=1,NY).
C (INPUT/OUTPUT)
C NX = THE SIZE OF THE ERROR( , ) GRID IN THE X
C DIRECTION. (INPUT)
C NY = THE SIZE OF THE ERROR( , ) GRID IN THE Y
C DIRECTION. (INPUT)
C CPNDFD(IX,JY) = THE NDFD MASK FROM THE MOS-2000 EXTERNAL
C RANDOM ACCESS FILE (IX=1,NY) (JY=1,NY) AT
C NOMINAL MESHLENGTH MESHE.
C 1 = WITHIN THE AREA; 0 = OUTSIDE. (INPUT)
C SEALND(IX,JY) = THE SEA/LAND (IX=1,NX) (JY=1,NY) AT
C NOMINAL MESHLENGTH MESHE.
C 9 = LAND; VALUES BELOW WATER. (INPUT)
C TELEV(IX,JY) = THE TERRAIN (IX=1,NX) (JY=1,NY) AT
C NOMINAL MESHLENGTH MESHE. (INPUT)
C NXE = X-EXTENT OF SEALND( ), AND CPNDFD( )
C AT MESH LENGTH MESHE. MUST = NX. (INPUT)
C NYE = Y-EXTENT OF SEALND( ), AND CPNDFD( )
C AT MESH LENGTH MESHE. MUST = NY. (INPUT)
C ITABLE(J,L,M) = HOLDS THE 2ND WORD IDS FOR UP TO JDIM
C PREDICTORS (J=1,JDIM) FOR UP TO IDIM
C EQUATIONS (M=1,IDIM) FOR LAND (L=1) AND
C WATER (L=2).. (INPUT)
C RTABLE(J,L,M) = HOLDS THE JDIM COEFFICIENTS (J=1,JDIM) AND
C THE CONSTANT (J=JDIM+1) FOR THE LAND (L=1)
C AND WATER (L=2) EQUATION FOR IDIM EQUATIONS
C (M=1,IDIM), EACH EQUATION PERTAINING TO A
C DIFFERENT VARIABLE. (INPUT)
C IDIM = THE MAXIMUM NUMBER OF PAIRS (LAND AND WATER)
C EQUATIONS. (INPUT)
C JDIM = THE MAXIMUM NUMBER OF TERMS IN THE EQUATIONS.
C (INPUT)
C M = THE NUMBER OF THE EQUATION. M WOULD VARY
C WITH WEATHER ELEMENT. (INPUT)
C RELVAR(J) = THE RADII OVER WHICH TO COMPUTE THE WEIGHTED
C VARIABILITY (J=1,4). (INPUT)
C VARTAB(J) = THE FACTOR TO INCREASE RELVAR( ) FOR LAND (J=1)
C AND WATER (J=2). (INPUT)
C ITERML = THE NUMBER OF TERMS IN THE LAND EQUATION
C EVALUATED. (INPUT/OUTPUT)
C ITERMW = THE NUMBER OF TERMS IN THE WATER EQUATION
C EVALUATED. (INPUT/OUTPUT)
C ISTOP = INCREMENTED BY 1 WHEN AN ERROR OCCURS.
C (INPUT/OUTPUT)
C IER = STATUS RETURN
C 0 = GOOD.
C 777 = CAN'T FIND A STATION WITHIN DISTANCE R.
C (OUTPUT)
C DISTSQ = DISTANCE (IN GRID UNITS) SQUARED BETWEEN
C A GRIDPOINT AND A STATION. (INTERNAL)
C R = THE THIRD PART OF THE ID WORD. THIS IS THE
C DISTANCE FROM THE GRIDPOINT TO DO THE SEARCH.
C (INTERNAL)
C VAR1 = THE VARIABILITY VARIABLE (SEE PURPOSE).
C (INTERNAL)
C JCOUNT = COUNTS THE NUMBER OF TIMES A RADIUS HAD TO BE
C INCREASED FOR THIS VARIABLE. (INTERNAL)
C TEMPX(J) = AN ARRAY TO HOLD ESTIMATES OF THE VALUE AT
C WITHHELD STATION K BASED ON THE VALUE AT EACH
C NEIGHBORING STATION, ADJUSTED BY THE NEIGHBOR'S
C LAPSE RATE AND VERTICAL SEPARATION BETWEEN
C THE TWO STATIONS (J=1,NSTA).
C (AUTOMATIC) (INTERNAL)
C TEMPW(J) = AN ARRAY TO HOLD QUADRATIC WEIGHTS COMPUTED
C FROM THE HORIZONTAL DISTANCES BETWEEN THE
C STATIONS AND THE RADIUS OF SEARCH (THE SAME
C WEIGHTING USED IN THE CORRECTION ALGORITHM)
C (J=1,NSTA). (AUTOMATIC) (INTERNAL)
C ICOUNT = COUNTS THE NUMBER OF GRIDPOINTS SET TO MISSING.
C (INTERNAL)
C RMULT = VALUE TO MUTIPLY BY THE RADIUS OF SEARCH FOR
C SPARSE DATA REGION (INPUT)
C 1 2 3 4 5 6 7 X
C
C NONSYSTEM SUBROUTINES CALLED
C NONE.
C
CHARACTER*8 CCALL(NSTA)
C
DIMENSION ID(4)
DIMENSION XP(NSTA),YP(NSTA),LTAG(NSTA),LNDSEA(NSTA),XDATA(NSTA),
1 XLAPSE(NSTA),ELEV(NSTA)
DIMENSION ERROR(NX,NY)
DIMENSION SEALND(NXE,NYE),CPNDFD(NXE,NYE),TELEV(NXE,NYE)
DIMENSION ITABLE(JDIM,2,IDIM),RTABLE(JDIM+1,2,IDIM)
DIMENSION RELVAR(4),VARTAB(2)
DIMENSION TEMPX(NSTA),TEMPW(NSTA)
C TEMPX( ) AND TEMPW( ) ARE AUTOMATIC ARRAYS LARGE ENOUGH TO
C HOLD THE MAXIMUM NUMBER OF POINTS SURROUNDING A STATION
C WITHIN ITS RADIUS OF SEARCH.
C
CALL TIMPR(KFILDO,KFILDO,'START VARILG ')
C
IER=0
IFIRST=0
ICOUNT=0
C ICOUNT COUNTS THE NUMBER OF POINTS SET TO MISSING.
IB=0
C IB = 1: INDICATES WHETHER A BACKUP VALUE OF RELVAR( ) HAS BEEN
C USED.
C IB = 2: INDICATES WHETHER R HAS BEEN INCREASED TO R*RMULT.
C
C THIS DEALS WITH BOTH WATER AND LAND EQUATIONS AND ALL TERMS
C IN EACH EQUATION. ERROR( , ) IS ACCUMULATED FOR ALL TERMS
C INVOLVING 0997, 1097, 1197, OR 1297 IN THE 2ND WORD.
C
DO 400 J=1,JDIM
C THERE ARE A MAXIMUM OF JDIM TERMS IN THE EQUATION.
DO 399 L=1,2
C L=1 FOR LAND, 2 = WATER.
C
IF(ITABLE(J,L,M)/10000.NE.1397.AND.
1 ITABLE(J,L,M)/10000.NE.1497.AND.
2 ITABLE(J,L,M)/10000.NE.1597.AND.
3 ITABLE(J,L,M)/10000.NE.1697)GO TO 399
C TRANSFER OUT IF THIS TERM IS NOT A VARIABILITY
C PREDICTOR.
C
C THIS IS A TERM TO EVALUATE.
C
IF(L.EQ.1)THEN
ITERML=ITERML+1
ELSE
ITERMW=ITERMW+1
ENDIF
C
C RETRIEVE R FROM THE ID.
C
R=ITABLE(J,L,M)-(ITABLE(J,L,M)/10000)*10000
C R IS THE DISTANCE TO DO THE SEARCH, TAKEN FROM THE ID.
RSQ=R*R+.01
C THE SMALL CONSTANT IS ADDED TO ASSURE A POINT IS NOT
C ELIMINATED BECAUSE OF ROUNDOFF, AND TO BE CONSISTENT WITH
C CLOS2
JCOUNT=0
C JCOUNT COUNTS THE NUMBER OF TIMES A RADIUS HAD TO BE
C INCREASED FOR THIS VARIABLE.
C
IF(NXE.NE.NX.OR.NYE.NE.NY)THEN
WRITE(KFILDO,101)NXE,NYE,NX,NY
101 FORMAT(/' ****NXE =',I5,' OR NYE =',I5,
1 ' DOES NOT MATCH NX =',I5,' OR NY =',I5,
2 ' AT 101 IN VARILG. FATAL ERROR.')
ISTOP=ISTOP+1
IER=777
GO TO 500
ENDIF
C
!$OMP PARALLEL DO DEFAULT(SHARED)
!$OMP& PRIVATE(IX,JY,K,DISTSQ,TEMPX,TEMPW,
!$OMP& VAR,NCOUNT,AVG,JJ,VAR1,ERR,MM,LL)
!$OMP& FIRSTPRIVATE(IB,R,RSQ)
C
DO 200 JY=1,NY
DO 199 IX=1,NX
C
IF(CPNDFD(IX,JY).LT..5)GO TO 199
C THIS GRIDPOINT IS OUTSIDE THE CLIPPING AREA.
C
C A POINT COULD HAVE BEEN SET TO MISSING.
C
IF(NINT(ERROR(IX,JY)).EQ.9999)GO TO 199
C
NCOUNT=0
VAR=0.
C
IF(L.EQ.1.AND.SEALND(IX,JY).GT.8.5)THEN
C
C THIS SECTION EVALUATES LAND POINTS.
C
120 DO 162 K=1,NSTA
IF(LNDSEA(K).LT.6)GO TO 162
C THIS IS A LAND GRIDPOINT, BUT A WATER STATION.
C ALLOWS THE INLAND WATER/LAND MIX TO BE USED.
C
IF(ABS(IX-XP(K)).GT.R)GO TO 162
IF(ABS(JY-YP(K)).GT.R)GO TO 162
C IN A LONG LIST OF STATIONS, THE ABOVE TWO TESTS SHOULD BE
C MORE EFFICIENT THAN ALWAYS CALCULATING THE DISTANCE.
C ALSO, THEY SHOULD RULE OUT MORE THAN THE TWO FOLLOWING
C TESTS.
IF(LTAG(K).NE.0)GO TO 162
C CHECKING LTAG(K) IS TANTAMOUNT TO CHECKING XDATA(K)
C FOR MISSING, AND ACCOUNTS FOR ANY TOSSED STATIONS ON
C THE LAST PASS.
DISTSQ=(IX-XP(K))**2+(JY-YP(K))**2
C
IF(DISTSQ.LT.RSQ)THEN
NCOUNT=NCOUNT+1
TEMPW(NCOUNT)=(RSQ-DISTSQ)/(RSQ+DISTSQ)
C WT CANNOT BE NEGATIVE, BECAUSE DISTSQ LT RSQ.
TEMPX(NCOUNT)=XDATA(K)+XLAPSE(K)*(TELEV(IX,JY)-ELEV(K))
C TEMPX(NCOUNT) IS THE VALUE IMPLIED BY STATION K AT
C THE GRIDPOINT.
VAR=VAR+TEMPX(NCOUNT)
C
CD IF(IX.EQ.300.AND.JY.EQ.400)THEN
CD WRITE(KFILDO,160)K,CCALL(K),IX,JY,DISTSQ,RSQ,XDATA(K),
CD 1 ELEV(K),TELEV(IX,JY),VAR,NCOUNT
CD160 FORMAT(/,' AT 160--K,CCALL(K),IX,JY,DISTSQ,RSQ,',
CD 1 'XDATA(K),ELEV(K),TELEV(IX,JY),VAR,NCOUNT',/,
CD 2 I6,1X,A8,2I6,2F8.1,4F8.2,I3)
CD ENDIF
C
ENDIF
162 CONTINUE
C
IF(NCOUNT.GT.1)THEN
C WHEN NCOUNT = 1, THERE WOULD BE ZERO VARIABILITY.
C THIS WOULD NOT GIVE A GOOD PICTURE, SO USE A BACKUP
C VALUE IN RELVAR( ).
AVG=VAR/NCOUNT
ELSE
C
IF(IFIRST.EQ.0)THEN
CD WRITE(KFILDO,163)R,IX,JY
CD163 FORMAT(/,' ####CANNOT FIND TWO STATIONS TO AVERAGE IN',
CD 1 ' VARILG WITHIN RADIUS R =',F5.1,' FOR',
CD 2 ' GRIDPOINT IX,JY',2I5,'. A LARGER RADIUS',
CD 3 ' WILL BE TRIED IN LAND.',/,
CD 4 ' THIS DIAGNOSTIC WILL NOT PRINT AGAIN.')
IFIRST=1
ENDIF
C
DO 1635 JJ=4,1,-1
C
IF(NINT(R).LE.NINT(RELVAR(JJ)*VARTAB(L)))THEN
C
IF(JJ.EQ.1)THEN
IB=2
C IB = 2 INDICATES R HAS BEEN INCREASED WITH RMULT.
R=R*RMULT
RSQ=R*R+.01
NCOUNT=0
VAR=0.
JCOUNT=JCOUNT+1
GO TO 120
ELSE
IB=1
C IB = 1 INDICATES R HAS BEEN INCREASED WITH RELVAR.
R=RELVAR(JJ-1)*VARTAB(L)
RSQ=R*R+.01
NCOUNT=0
VAR=0.
JCOUNT=JCOUNT+1
CD WRITE(KFILDO,1634)IX,JY,R,RSQ
CD1634 FORMAT(/,' AT 1634 IN VARILG--IX,JY,R,RSQ',2I5,2F10.2)
GO TO 120
ENDIF
C
ELSE
IF(IB.EQ.2)THEN
WRITE(KFILDO,1633)IX,JY
1633 FORMAT(' *****LAND POINT IX, JY =',2I6,' CANNOT BE',
1 ' EVALUATED AND SET TO MISSING IN VARILG.',
2 ' MAY MEAN MISSING DATA',
3 ' IN THE ANALYSIS.')
ERROR(IX,JY)=9999.
ICOUNT=ICOUNT+1
GO TO 199
ENDIF
ENDIF
C
1635 CONTINUE
C
ENDIF
C
IF(IB.NE.0)THEN
C A LARGER VALUE OF R HAS BEEN USED THAN THE ID SPECIFIES.
C THEREFORE, RESET R AND RSQ TO THE CORRECT VALUES.
R=ITABLE(J,L,M)-(ITABLE(J,L,M)/10000)*10000
C R IS THE DISTANCE TO DO THE SEARCH, TAKEN FROM THE ID.
RSQ=R*R+.01
C THE SMALL CONSTANT IS ADDED TO ASSURE A POINT IS NOT
C ELIMINATED BECAUSE OF ROUNDOFF, AND TO BE CONSISTENT WITH
C CLOS2
IB=0
CD WRITE(KFILDO,1636)IX,JY,R,RSQ
CD1636 FORMAT(' AT 1636 IN VARILG--IX,JY,R,RSQ',2I5,2F10.2)
ENDIF
C
VAR1=0.
C
DO 164 LL=1,NCOUNT
VAR1=VAR1+ABS(TEMPX(LL)-AVG)*TEMPW(LL)
C
CD IF(IX.EQ.300.AND.JY.EQ.400)THEN
CD WRITE(KFILDO,1637)J,L,M,IX,JY,VAR1,TEMPX(LL),TEMPW(LL),AVG
CD1637 FORMAT(/,' AT 1637 IN VARILG--J,L,M,IX,JY,VAR1,',
CD 1 'TEMPX(LL),TEMPW(LL),AVG',/,5I6,4F10.4)
CD ENDIF
C
164 CONTINUE
C
ERR=VAR1/NCOUNT
C NCOUNT HAS BEEN PREVIOUSLY CHECKED FOR GT 1.
ERROR(IX,JY)=ERROR(IX,JY)+ERR*RTABLE(J,L,M)
C
CCC IF(IX.EQ.300.AND.JY.EQ.400)THEN
CCC WRITE(KFILDO,165)J,L,M,IX,JY,ERR,RTABLE(J,L,M),
CCC 1 ERROR(IX,JY)
CCC 165 FORMAT(/,' AT 165 IN VARILG--J,L,M,IX,JY,ERR,',
CCC 1 'RTABLE(J,L,M),ERROR(IX,JY)',/,
CCC 2 5I6,3F10.4)
CCC ENDIF
C
ELSEIF(L.EQ.2.AND.SEALND(IX,JY).LT.8.5)THEN
C
C THIS SECTION EVALUATES WATER POINTS.
C
170 DO 172 K=1,NSTA
IF(LNDSEA(K).GT.6)GO TO 172
C THIS IS A WATER GRIDPOINT, BUT A LAND STATION.
C ALLOWS THE INLAND WATER/LAND MIX TO BE USED.
C
IF(ABS(IX-XP(K)).GT.R)GO TO 172
IF(ABS(JY-YP(K)).GT.R)GO TO 172
C IN A LONG LIST OF STATIONS, THE ABOVE TWO TESTS SHOULD BE
C MORE EFFICIENT THAN ALWAYS CALCULATING THE DISTANCE.
C ALSO, THEY SHOULD RULE OUT MORE THAN THE TWO FOLLOWING
C TESTS.
IF(LTAG(K).NE.0)GO TO 172
C CHECKING LTAG(K) IS TANTAMOUNT TO CHECKING XDATA(K)
C FOR MISSING, AND ACCOUNTS FOR ANY TOSSED STATIONS ON
C THE LAST PASS.
DISTSQ=(IX-XP(K))**2+(JY-YP(K))**2
C
IF(DISTSQ.LT.RSQ)THEN
NCOUNT=NCOUNT+1
TEMPW(NCOUNT)=(RSQ-DISTSQ)/(RSQ+DISTSQ)
C WT CANNOT BE NEGATIVE, BECAUSE DISTSQ LT RSQ.
TEMPX(NCOUNT)=XDATA(K)+XLAPSE(K)*(TELEV(IX,JY)-ELEV(K))
C TEMPX(NCOUNT) IS THE VALUE IMPLIED BY STATION K AT
C THE GRIDPOINT.
VAR=VAR+TEMPX(NCOUNT)
C
ENDIF
172 CONTINUE
C
IF(NCOUNT.GT.1)THEN
C WHEN NCOUNT = 1, THERE WILL BE ZERO VARIABILITY.
C THIS DOES NOT GIVE A GOOD PICTURE, SO SET TO 9999.
AVG=VAR/NCOUNT
ELSE
C
IF(IFIRST.EQ.0)THEN
CD WRITE(KFILDO,173)R,IX,JY
CD173 FORMAT(/,' ####CANNOT FIND TWO STATIONS TO AVERAGE IN',
CD 1 ' VARILG WITHIN RADIUS R =',F5.1,' FOR',
CD 2 ' GRIDPOINT IX,JY',2I5,'. A LARGER RADIUS',
CD 3 ' WILL BE TRIED IN WATER.',/,
CD 4 ' THIS DIAGNOSTIC WILL NOT PRINT AGAIN.')
IFIRST=1
ENDIF
C
DO 1735 JJ=4,1,-1
C
IF(NINT(R).LE.NINT(RELVAR(JJ)*VARTAB(L)))THEN
C
IF(JJ.EQ.1)THEN
IB=2
C IB = 2 INDICATES R HAS BEEN INCREASED WITH RMULT.
R=R*RMULT
RSQ=R*R+.01
NCOUNT=0
VAR=0.
JCOUNT=JCOUNT+1
GO TO 170
ELSE
IB=1
C IB = 1 INDICATES R HAS BEEN INCREASED WITH RELVAR.
R=RELVAR(JJ-1)*VARTAB(L)
RSQ=R*R+.01
NCOUNT=0
VAR=0.
JCOUNT=JCOUNT+1
CD WRITE(KFILDO,1734)IX,JY,R,RSQ
CD1734 FORMAT(/,' AT 1734 IN VARILG--IX,JY,R,RSQ',
CD 1 2I5,2F10.2)
GO TO 170
ENDIF
C
ELSE
IF(IB.EQ.2)THEN
WRITE(KFILDO,1733)IX,JY
1733 FORMAT(' *****WATER POINT IX, JY =',2I6,' CANNOT BE',
1 ' EVALUATED AND SET TO MISSING IN VARILG.',
2 ' MAY MEAN MISSING DATA',
3 ' IN THE ANALYSIS.')
ERROR(IX,JY)=9999.
ICOUNT=ICOUNT+1
GO TO 199
ENDIF
ENDIF
C
1735 CONTINUE
C
ENDIF
C
IF(IB.NE.0)THEN
C A LARGER VALUE OF R HAS BEEN USED THAN THE ID SPECIFIES.
C THEREFORE, RESET R AND RSQ TO THE CORRECT VALUES.
R=ITABLE(J,L,M)-(ITABLE(J,L,M)/10000)*10000
C R IS THE DISTANCE TO DO THE SEARCH, TAKEN FROM THE ID.
RSQ=R*R+.01
C THE SMALL CONSTANT IS ADDED TO ASSURE A POINT IS NOT
C ELIMINATED BECAUSE OF ROUNDOFF, AND TO BE CONSISTENT WITH
C CLOS2
IB=0
CD WRITE(KFILDO,1737)IX,JY,R,RSQ
CD1737 FORMAT(' AT 1737 IN VARIG--IX,JY,R,RSQ',2I5,2F10.2)
ENDIF
C
VAR1=0.
C
DO 182 LL=1,NCOUNT
VAR1=VAR1+ABS(TEMPX(LL)-AVG)*TEMPW(LL)
182 CONTINUE
C
ERR=VAR1/NCOUNT
C NCOUNT HAS BEEN PREVIOUSLY CHECKED FOR GT 1.
ERROR(IX,JY)=ERROR(IX,JY)+ERR*RTABLE(J,L,M)
C
IF(IX.EQ.300.AND.JY.EQ.400)THEN
WRITE(KFILDO,185)J,L,M,IX,JY,ERR,RTABLE(J,L,M),
1 ERROR(IX,JY)
185 FORMAT(/,' AT 185 IN VARILG--J,L,M,IX,JY,ERR,',
1 'RTABLE(J,L,M),ERROR(IX,JY)',/,
2 5I6,3F10.4)
ENDIF
C
ENDIF
C
199 CONTINUE
200 CONTINUE
C
!$OMP END PARALLEL DO
C
IF(L.EQ.1)THEN
C
IF(JCOUNT.EQ.0)THEN
WRITE(KFILDO,250)ID(1),ITABLE(J,L,M)
250 FORMAT('THE RADIUS OF SEARCH DID NOT HAVE TO BE',
1 ' INCREASED IN VARILG FOR THIS VARIABLE ',2I10,
2 ' FOR LAND POINTS.')
ELSE
WRITE(KFILDO,251)JCOUNT,ID(1),ITABLE(J,L,M)
251 FORMAT('THE RADIUS OF SEARCH HAD TO BE INCREASED',
1 I8,' TIMES IN VARILG FOR THIS VARIABLE ',2I10,
2 ' FOR LAND POINTS.')
ENDIF
C
ELSE
IF(JCOUNT.EQ.0)THEN
WRITE(KFILDO,260)ID(1),ITABLE(J,L,M)
260 FORMAT('THE RADIUS OF SEARCH DID NOT HAVE TO BE',
1 ' INCREASED IN VARILG FOR THIS VARIABLE ',2I10,
2 ' FOR WATER POINTS.')
ELSE
WRITE(KFILDO,261)JCOUNT,ID(1),ITABLE(J,L,M)
261 FORMAT('THE RADIUS OF SEARCH HAD TO BE INCREASED',
1 I8,' TIMES IN VARILG FOR THIS VARIABLE ',2I10,
2 ' FOR WATER POINTS.')
ENDIF
C
ENDIF
C
399 CONTINUE
400 CONTINUE
C
IF(ICOUNT.GT.0)THEN
WRITE(KFILDO,410)ICOUNT
410 FORMAT(/' ####NUMBER OF POINTS SET TO MISSING IN VARILG =',I9)
ENDIF
C
CALL TIMPR(KFILDO,KFILDO,'END VARILG ')
C
500 RETURN
END