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