SUBROUTINE WTHOL2(KFILDO,IP24,ID,NDATE,CCALL,NAME,XDATA,ND1, 1 LTAG,LNDSEA,ELEV,WHOLD,LTAGWH,LOCWH,XP,YP, 2 DUMCAL,XLAPSE,NSTA, 3 CPNDFD,SEALND,TELEV,NXE,NYE,MESHE,NCLIPY, 4 NX,NY,MESH,SEED,NWITH,IEXTRA, 5 NROUGH,RELVAR,VARTAB, 6 ERRANL,ID2,IDIM,L3264B,L3264W,ISTOP,IER) C C NOVEMBER 2007 GLAHN MDL MOS-2000 C MODIFIED FROM WTHOL1 C NOVEMBER 2007 GLAHN 2ND EDITION C NOVEMBER 2007 GLAHN REMOVED DIMENSION ON RCLOS2 C NOVEMBER 2007 GLAHN RECOMPUTED X WHEN Y WAS NOT C USABLE; REQUIRED GRIDPOINT BE LAND C DECEMBER 2007 GLAHN MODIFIED FOR LIO AND DIFFERENT C DATA TYPES; CHANGED VALUES OF C RVARI( ) AND RELVAR( ) C DECEMBER 2007 GLAHN MINIMIZED MISSINGS BY USING NEXT C LARGER RADIUS FOR VARI AND ELVAR; C EXCLUDED CANADIAN STATIONS; ADDED C NAME( ) C DECEMBER 2007 GLAHN RCLOS2 SET 45. VICE 35. C DECEMBER 2007 GLAHN ADDED ID2, IDEM, JDEM; MADE RADIUS C 27.5 A WHOLE NUMBER 27.0 C DECEMBER 2007 GLAHN ADDED ELD, ELDS, ELDXCL TO ERRANL C DECEMBER 2007 GLAHN ADDED DEFINITION OF ID2(20) (12/29) C DECEMBER 2007 GLAHN ADDED NROUGH,RCLOS2 RVARI,RDENS, C RELVAR TO CALL C JANUARY 2008 GLAHN ELIMINATED DENS, ELVAR C FEBRUARY 2008 GLAHN AFTER CALL TO VARI AND VARIW, C TRANSFERRED TO 130 ONLY FOR M=1 C FEBRUARY 2008 GLAHN ELIMINATED RDENS C FEBRUARY 2008 GLAHN NAME OF VARIW CHANGED TO VARIL C ADDED VARIW FOR VARIABLES 3 AND 4 C DECEMBER 2008 GLAHN REMOVED $NOSTANDARD INTRINSICS C BECAUSE OF IBM COMPILER C DECEMBER 2008 GLAHN HP FUNCTION RAN REPLACED WITH IBM C 'RANDOM_NUMBER' ADN 'RANDOM_SEED' C JANUARY 2009 GLAHN RANDOM NUMBER GENERATOR NOW RAND( ) C FEBRUARY 2009 GLAHN CALLED SRAND(SEED) ONLY WHEN JFIRST =0 C FEBRUARY 2009 GLAHN MINOR COMMENT CHANGES C MARCH 2009 GLAHN REMOVED IER CHECKING CODE FOLLOWING C CALLS TO VARI, VARIW, AND VARIL; C ELIMINATED RCLOS2 FROM CALL; ELIMINATED C RVARI ALTOGETHER; ADDED VARTAB( ) C AUGUST 2009 GLAHN COMMA PUT IN FORMAT 901 C DECEMBER 2009 GLAHN IMPROVED COMMENT FOR ERRANL( ,18) C JUNE 2014 GLAHN SET IER = 0;INSERTED C CALL W3TAGE('155') BEFORE STOPS C DECEMBER 2014 GHIRARDELLI MODIFIED OPEN FOR OPERATIONS C C C PURPOSE C TO WITHHOLD STATIONS FROM THE DATA TO ANALYZE FOR U155. C THE DATA WITHHELD ARE SAVED IN WHOLD( ). LTAG( ) IS C MODIFIED ACCORDING TO THE WITHHELD DATA, AND THE VALUES C SAVED IN LTAGWH( ) FOR RESTORATION. C SELECTION FOR WITHHOLDING IS TO RANDOMLY PICK A C POSITION ON THE GRID IN THE X DIRECTION AND ALSO IN C THE Y DIRECTION AND FIND THE CLOSEST STATION WITH C NON-MISSING DATA. ANY GRIDPOINT OUTSIDE THE NDFD C GRID IS NOT USED; THIS IS DONE USING CPNDFD( ). C MESH MUST EQUAL MESHE FOR CLIPPING TO BE DONE. C C NOTE (1): THIS ROUTINE IS FOR ASSESSING THE ACCURACY OF C THE ANALYSIS. NWITH STATIONS ARE WITHHELD FROM C THE ANALYSIS, AND ANOTHER IEXTRA USED IN THE C ANALYSIS ARE SAVED. THESE DATA ARE WRITTEN FOR C LATER ANALYSIS. C NOTE (2): THIS ROUTINE WILL WORK FOR LAND OR WATER, C POINTS, DEPENDING ON THE VARIABLE LIO SET C BY DATA STATEMENT TO 9 OR 0, RESPECTIVELY. C NOTE (3): THE ORIGINAL SELECTION OF A STATION CLOSE C TO A RANDOM X/Y POINT IS BY ONE OF TWO TYPES-- C LIO = 9, USE LNDSEA( ) = 6 AND 9 C LIO = 0, USE LNDSEA( ) = 0 AND 3. C NOTE (4): NO DISTINCTION IS MADE BETWEEN OCEAN AND C INLAND WATER BECAUSE THE POINTS ARE TOO FEW C TO HAVE SEPARATE RELATIONSHIPS. C NOTE (5): ONLY GRIDPOINTS WITHIN THE NDFD CLIPPING C GRID CPNDFD( ) ARE USED. C NOTE (6): CANADIAN STATIONS WITH CALL LETTERS STARTING C WITH 'CW' AND 'CY' ARE NOT USED. C C DATA SET USE C KFILDO - UNIT NUMBER OF OUTPUT (PRINT) FILE. (OUTPUT) C IP24 - UNIT NUMBER TO WRITE LIST OF WITHHELD STATIONS. C (OUTPUT) C C VARIABLES C KFILDO = UNIT NUMBER OF OUTPUT (PRINT) FILE. (INPUT) C IP24 = UNIT NUMBER FOR WRITING LIST OF WITHHELD C STATIONS AND DATA. (INPUT) C ID(J) = THE VARIABLE ID'S BEING DEALT WITH(J=1,4). C (INPUT) C NDATE = THE DATE/TIME OF THE RUN. (INPUT) C CCALL(K) = CALL LETTERS OF STATIONS (K=1,NSTA). C (CHARACTER*8) (INPUT) C NAME(K) = NAMES OF STATIONS (K=1,NSTA). (CHARACTER*20) C (INPUT) C XDATA(K) = ON INPUT, HOLDS THE FULL SET OF DATA C (K=1,NSTA). (INPUT) C ND1 = FIRST DIMENSION OF XDATA( ). (INPUT) C LTAG(K) = DENOTES USE OF DATA IN DATA(K) FOR STATION K C (K=1,NSTA). C 0 = USE DATA. C 1 = STATION OUTSIDE RADIUS OF INFLUENCE FOR C AREA BEING ANALYZED OR MISSING DATUM. C 2 = STATION LOCATION UNKNOWN. C (INPUT/OUTPUT) 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 ELEV(K) = ELEVATION OF STATIONS (K=1,NSTA). (INPUT) C WHOLD(K) = ARRAY TO HOLD DATA WITHHELD (K=1,NSTA). C ON EXIT, A VALUE WITHHELD WILL BE IN WHOLD( ); C OTHERWISE, WHOLD( ) = -9999. (OUTPUT). C LTAGWH(K) = THE VALUES OF LTAG( ) OF THE WITHHELD STATIONS C ARE SAVED SO THEY CAN BE RESTORED (K=1,NSTA). C (OUTPUT) C LOCWH(K) = THE LOCATIONS OF THE WITHHELD STATIONS IN C THE LIST IN THE ORDER OF WITHHOLDING. C (OUTPUT) C XP(K) = THE X POSITION FOR STATION K (K=1,NSTA) ON C THE ANALYSIS GRID AREA AT THE CURRENT NOMINAL C GRID MESH LENGTH MESH. (INPUT) C YP(K) = THE Y POSITION FOR STATION K (K=1,NSTA) ON C THE ANALYSIS GRID AREA AT THE CURRENT NOMINAL C GRID MESH LENGTH MESH. (INPUT) C DUMCAL(K) = DUMMY CALL LETTERS OF THE STATIONS PROCESSED C (K=1,NWITH+IEXTRA). (CHARACTER*8) (OUTPUT) C (AUTOMATIC) 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 OR LOCATIONS BEING DEALT C WITH. (INPUT) C CPNDFD(J) = THE NDFD MASK FROM THE MOS-2000 EXTERNAL C RANDOM ACCESS FILE (J=1,NXE*NYE) AT NOMINAL C MESHLENGTH MESHE. (INPUT) C SEALND(J) = THE LAND/SEA MASK (J=1,NXE*NYE). (INPUT) C 0 = OCEAN WATER GRIDPOINTS; C 3 = INLAND WATER GRIDPOINTS. C 9 = LAND GRIDPOINTS. C TELEV(J) = THE TERRAIN ELEVATION (J=1,NXE*NYE). (INPUT) C NXE = X-EXTENT OF TELEV( ), SEALND( ), AND CPNDFD( ) C AT MESH LENGTH MESHE. (INPUT) C NYE = Y-EXTENT OF TELEV( ), SEALND( ), AND CPNDFD( ) C AT MESH LENGTH MESHE. (INPUT) C MESHE = THE NOMINAL MESH LENGTH OF THE TERRAIN GRID. C IT IS MANDATORY THE GRID AVAILABLE IS OF THIS C MESH SIZE AND COVER THE SAME AREA SPECIFIED C BY NXL BY NYL, EVEN IF MESHE IS NOT EQUAL C TO MESHB. (INPUT) C NCLIPY = 1 WHEN THE NDFD MASK GRID IS AVAILABLE AND C IN CPNDFD( ). C 0 OTHERWISE. C (INPUT) C NX = THE SIZE OF THE GRID IN THE X DIRECTION. C (INPUT) C NY = THE SIZE OF THE GRID IN THE Y DIRECTION. C (INPUT) C MESH = THE MESH LENGTH OF THE GRID TO WHICH THE C POSITIONS XP( ) AND YP( ) REFER. (INPUT) C SEED = SEED FOR THE RANDOM NUMBER GENERATOR. C (INPUT/OUTPUT) C NWITH = NUMBER OF STATIONS TO WITHHOLD. (INPUT) C IEXTRA = THE NUMBER OF STATIONS TO PROCESS BUT TO C NOT WITHHOLD. (INPUT) C NROUGH(J) = THE RADII OVER WHICH TO COMPUTE THE ROUGHNESS C (J=1,4). SET BY DATA STATEMENT IIN U405A. C THE VALUES CURRENTLY ARE IN ORDER 8, 4, 2, 1. C (INPUT) C RELVAR(J) = THE RADII OVER WHICH TO COMPUTE THE VARIABLILTY C VARIABLES IN TERMS OF GRIDLENGTHS (J=1,4). C ALSO, RELVAR(1) IS THE DISTANCE OVER WHICH TO C SEARCH FOR THE CLOSEST TWO STATIONS. (INPUT) C VARTAB(J) = THE MULTIPLYING FACTORS FOR RELVAR( ) FOR C LAND(J=1) AND WATER(2) FOR AREA BEING DEALT C WITH. (INPUT) C ERRANL(K,J) = THE VALUES FOR THE ERROR ANALYSIS C (K=1,NWITH+IEXTRA) (J=1,20). THE COLUMNS ARE: C 1,2 -- THE DISTANCE TO THE TWO CLOSEST STATIONS C 3,4 -- DATA VARIABILITY, WITH LAPSE BUT NOT WEIGHTS C 5-8 -- THE 4 TERRAIN ROUGHNESS VARIABLES C 9-12 -- DATA VARIABILITY, NO LAPSE OR WEIGHTS C 13-16 -- DATA VARIABILITY, WITH LAPSE AND WEIGHTS C AND ITS NEIGHBORS, ACCOUNTING FOR XLAPSE( ) C 17 -- ABSOLUTE DIFFERENCE IN ELEVATION BETWEEN C STATION AND ITS CLOSEST NEIGHBOR C 18 -- ABSOLUTE DIFFERENCE BETWEEN STATION DATA C VALUE AND VALUE OF THE CLOSEST NEIGHBOR C AFTER LAPSE AT THE CLOSEST STATION IS C APPLIED C 19 -- THE PRODUCT OF 17 AND 1 C 20 -- THE "ERROR," THE DIFFERENCE BETWEEN THE C ACTUAL VALUE AND THE INTERPOLATED VALUE C FROM THE ANALYSIS C (OUTPUT) C ID2(M) = THE 2ND WORD ID OF THE VARIABLE M USED IN C ERRANL( , ) FOR THE ERROR DATA (M=1,IDIM). C IT IS COMPOSED OF 3 PARTS XXX97YYYY: C (1) XXX = NUMBER IN SEQUENCE 1-20 (OR MORE), C (2) CONSTANT = 97, C (3) YYYY = THE VALUE USED IN COMPUTING THE C VARIABLE (E.G.,35 WHEN THE VARIABILITY C HAS BEEN COMPUTED OVER 35 GRIDLENGTHS C FROM THE STATION OR GRIDPOINT). C (OUTPUT) C IDIM = THE NUMBER OF VARIABLES IN ID2( ) AND C ERRANL( , ) = IDEM. (INPUT) C L3264B = INTEGER WORD LENGTH IN BITS OF MACHINE BEING C USED (EITHER 32 OR 64). (INPUT). C L3264W = NUMBER OF WORDS IN 64 BITS (EITHER 1 OR 2). C (INPUT) C ISTOP = INCREMENTED BY 1 WHEN AN ERROR OCCURS. C (INPUT/OUTPUT) C IER = 0 = GOOD RETURN. (OUTPUT) C IFIRST = 0 ON FIRST ENTRY. IF A CLIPPING GRID IS NOT C AVAILABLE, A DIAGNOSTIC IS WRITTEN AND IFIRST C IS SET OT 1. IT IS SAVED FROM ENTRY TO ENTRY. C THE DIAGNOSTIC SHOULD BE WRITTEN ONLY ONCE C PER RUN. (INTERNAL) (SAVED) C JFIRST = 0 ON FIRST ENTRY. AFTER RAN IS CALLED 10 TIMES C IT IS SET TO 1. (INTERNAL) (SAVED) C MFIRST = 0 ON FIRST ENTRY. AFTER IP24 IS DELETED AND C REOPENED, IT IS SET TO 1. (INTERNAL) (SAVED) C WSTA(K) = CALL LETTERS OF THE STATIONS PROCESSED C (K=1,NWITH+IEXTRA). (CHARACTER*8) (INTERNAL) C (AUTOMATIC) C LIO = THE VARIABLE SET BY DATA STATEMENT TO SPECIFY C WHETHER THIS RUN IS FOR LAND (=9), OR WATER C (=0) POINTS. (INTERNAL) C LIOVAR = SET TO 1 FOR LAND STATIONS, TO 2 FOR WATER. C FOR ACCESS ITO VARTAB( ). (INTERNAL) C RCLOS2 = THE SEARCH RADIUS TO FURNISH TO CLOS2, WHERE C RCLOS2=RELVAR(1)*LIOVAR). (INTERNAL) C 1 2 3 4 5 6 7 X C C NONSYSTEM SUBROUTINES CALLED C ROUGH, CLOS2, VARI, VARIL C CHARACTER*8 FLNAM CHARACTER*8 CCALL(NSTA) CHARACTER*8 WSTA(NWITH+IEXTRA) C WSTA( ) IS AN AUTOMATIC ARRAY. CHARACTER*8 DUMCAL(NWITH+IEXTRA) CHARACTER*20 NAME(ND1) C DIMENSION ID(4) DIMENSION XDATA(NSTA),LTAG(NSTA),LNDSEA(NSTA),WHOLD(NSTA), 1 LTAGWH(NSTA),XP(NSTA),YP(NSTA),ELEV(NSTA),XLAPSE(NSTA) DIMENSION SEALND(NXE,NYE),TELEV(NXE,NYE),CPNDFD(NXE,NYE) DIMENSION LOCWH(NWITH+IEXTRA) DIMENSION ERRANL(NWITH+IEXTRA,IDIM),ID2(IDIM) DIMENSION NROUGH(4),RELVAR(4),VARTAB(2) C DATA IFIRST/0/, 1 JFIRST/0/, 2 MFIRST/0/ DATA LIO/9/ C****************************************************************** C FOR LIO = 9, THIS RUN IS FOR LAND POINTS. CHANGE IT TO 0 C FOR WATER POINTS. C****************************************************************** C SAVE IFIRST,JFIRST,MFIRST C IER=0 C WRITE(KFILDO,105)NROUGH,RELVAR,VARTAB 105 FORMAT(' AT 105 IN WTHOL2--NROUGH,RELVAR,VARTAB',/, 1 4I4,6F10.1) C C CALL RAN 10 TIMES THE FIRST CALL TO GET PAST ANY PROBLEMS WITH C USING AN INAPPROPRIATE INITIAL SEED . C CALL TIMPR(KFILDO,KFILDO,'START WTHOL2 ') C D WRITE(KFILDO,110)NWITH,(ID(J),J=1,4),NDATE D110 FORMAT(/,I4,' STATIONS ARE BEING WITHHELD PER ANALYSIS FOR', D 1 ' VARIABLE',3I10.9,I10.3,' FOR DATE',I12) C C CALL RAN 10 TIMES TO GET PAST ANY PROBLEMS WITH USING AN C INAPPROPRIATE INITIAL SEED THE FIRST TIME ENTERED. C IF(JFIRST.EQ.0)THEN C CALL SRAND(SEED) C THIS SETS THE SEED FOR THE RANDOM NUMBER GENERATOR. C DO 120 J=1,10 X=RAND() 120 CONTINUE C JFIRST=1 ENDIF C C INITIALIZE WHOLD( ) TO -9999. C DO 125 K=1,NSTA WHOLD(K)=-9999. 125 CONTINUE C IF(LIO.EQ.9)LIOVAR=1 IF(LIO.EQ.0)LIOVAR=2 RCLOS2=RELVAR(1)*VARTAB(LIOVAR) C RCLOS2 TAKES INTO ACCOUNT LAND OR WATER. RELVAR(1) IS C USED SO THAT IF A CLOSEST STATION IS FOUND, THEN A C CLOSEST STATION CAN BE FOUND IN OTHER ROUTINES. C C WITHHOLD NWITH STATIONS AND COMPUTE VARIABLES FOR C IEXTRA OTHERS. C DO 180 J=1,NWITH+IEXTRA C C PICK THE X POSITION. c 130 X=RAND() X=X*NX C X IS THE LOCATION IN THE GRID IN THE X DIRECTION. C D WRITE(KFILDO,132)J,SEED,X D132 FORMAT(' AT 132 IN WTHOL2--J,SEED,X',I4,F20.9,F20.9) C IF(X.LT.1..OR.X.GE.NX)GO TO 130 C THIS IS FOR SAFETY. ZERO COULD OCCUR. C C PICK THE Y POSITION. C Y=RAND() Y=Y*NY C Y IS THE LOCATION IN THE GRID IN THE Y DIRECTION. C D WRITE(KFILDO,142)J,SEED,Y D142 FORMAT(' AT 142 IN WTHOL2--J,SEED,Y',I4,F20.9,F20.9) C IF(Y.LT.1..OR.Y.GE.NY)GO TO 130 C THIS IS FOR SAFETY. ZERO COULD OCCUR. C RECALCULATE X AS WELL AS Y OR THIS CAN BIAS RESULTS. C CHANGED 11/28/07 C C THROW OUT A GRIDPOINT THAT IS NOT WITHIN THE NDFD GRID C AND IS NOT THE TYPE SPECIFIED BY LIO. C IF(NCLIPY.EQ.0)THEN C IF(IFIRST.EQ.0)THEN WRITE(KFILDO,145) 145 FORMAT(/,' ****THE CLIPPING GRID CPNDFD( ) IS NOT', 1 ' AVAILABLE. THIS WILL LIKELY BIAS RESULTS.', 2 ' PROCEEDING AT 145 IN WTHOL2.') IFIRST=1 ISTOP=ISTOP+1 ENDIF C ELSEIF(MESH.NE.MESHE)THEN C WRITE(KFILDO,146) 146 FORMAT(/,' ****THE CLIPPING GRID CPNDFD( ) IS NOT', 1 ' AT THE SAME MESH LENGTH AS THE GRID BEING', 2 ' USED. CLIPPING NOT DONE.',/, 3 ' FATAL ERROR. STOP IN WTHOL2', 2 'AT 146.') CALL W3TAGE('WTHOL2') STOP 146 C ELSE D WRITE(KFILDO,147)X,Y,CPNDFD(NINT(X),NINT(Y)) D147 FORMAT(' IN WTHOL2 AT 147--X,Y,CPNDFD(NINT(X),NINT(Y))', D 1 3F10.3) C IF(CPNDFD(NINT(X),NINT(Y)).EQ.0.)GO TO 130 C THE CLOSEST GRIDPOINT TO THE RANDOM POINT SELECTED C IS NOT WITHIN THE NDFD AREA. C C MAKE SURE TEH GRIDPOINT AND THE STATION ARE OF THE C SAME TYPE. C IF(LIO.EQ.9)THEN IF(SEALND(NINT(X),NINT(Y)).LT.6.5)GO TO 130 ELSEIF(LIO.EQ.0)THEN IF(SEALND(NINT(X),NINT(Y)).GT.3.5)GO TO 130 ELSE WRITE(KFILDO,1475)LIO 1475 FORMAT(/,' ****LIO =',I4,' NOT ONE OF THE VALUES TWO', 1 '0 OR 9. FATAL ERROR. STOP IN WTHOL2', 2 'AT 1475.') CALL W3TAGE('WTHOL2') STOP 1475 ENDIF C ENDIF C C FIND THE CLOSEST STATION OF THE SAME TYPE WITH NON-MISSING C DATA FOR THE DATA TO BE ANALYZED. C DISTSQ=99999999. KSTA=9999 C C DO 150 K=1,NSTA C CCCD WRITE(KFILDO,1477)K,LNDSEA(K),LTAG(K),XDATA(K) CCCD1477 FORMAT(/'AT 1477 IN WTHOL2--K,LNDSEA(K),LTAG(K),XDATA(K)', CCCD 1 3I6,F10.1) C IF(XDATA(K).EQ.9999.)GO TO 150 C DON'T DEAL WITH A MISSING DATUM. C IF(LTAG(K).NE.0)GO TO 150 C LTAG( ) = 0 MEANS THIS POINT IS NOT USED. C IF(LIO.EQ.9.AND.LNDSEA(K).LT.6)GO TO 150 IF(LIO.EQ.0.AND.LNDSEA(K).GT.3)GO TO 150 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). C IF(XP(K).LT.1..OR.XP(K).GT.NX)GO TO 150 IF(YP(K).LT.1..OR.YP(K).GT.NY)GO TO 150 C DON'T WITHHOLD STATIONS OUTSIDE GRID. THIS INSURES C 4 GRIDPOINTS ARE AVAILABLE FOR LINEAR INTERPOLATION. C CCCD WRITE(KFILDO,1499)K,CCALL(K),NAME(K), CCCD 1 KSTA,CCALL(KSTA),DISTSQ CCCD1499 FORMAT(' AT 1499 IN WTHOL2--K,CCALL(K),NAME(K),', CCCD 1 'KSTA,CCALL(KSTA),DISTSQ', CCCD 2 I6,2X,A8,2X,A20,i6,2x,a8,F8.1) C DQ=(X-XP(K))**2+(Y-YP(K))**2 C IF(DQ.LT.DISTSQ)THEN DISTSQ=DQ KSTA=K ENDIF C 150 CONTINUE C IF(KSTA.EQ.9999)THEN WRITE(KFILDO,151)X,Y 151 FORMAT(/,' ****COULD NOT FIND STATION OF CORRECT TYPE FOR', 1 ' RANDOM X,Y LOCATION =',2F7.2, 2 '. THIS SHOULD NOT HAPPEN. TRY ANOTHER POINT.') ISTOP=ISTOP+1 GO TO 130 ELSE C C DO NOT WITHHOLD CANADIAN STATIONS. THESE ARE STATIONS C WITH CALL LETTERS STARTING WITH 'CW' OR 'CY'. THERE C MAY BE ONE OR TOW US STATIONS THAT MEET THAT DEFINITION, C BUT THEY ARE VERY RARE. C IF(CCALL(KSTA)(1:2).EQ.'CW'.OR. 1 CCALL(KSTA)(1:2).EQ.'CY')THEN WRITE(KFILDO,1515)CCALL(KSTA),NAME(KSTA) 1515 FORMAT(/' ####CANADIAN STATION ',A8,A20, 1 ' NOT BEING WITHHELD. TRY ANOTHER POINT.') GO TO 130 ENDIF C C THIS IS A POINT TO PROCESS. C DIST=SQRT(DISTSQ) LOCWH(J)=KSTA C D WRITE(KFILDO,152)CCALL(KSTA),DIST D152 FORMAT(/,' AT 152 IN WTHOL2--CCALL(KSTA),DIST ',A8,F8.2) C C FIND THE TWO CLOSEST STATIONS WITHIN THE SEARCH RADIUS C RELVAR(1)*VARTAB(. C CALL CLOS2(KFILDO,CCALL,XP,YP,LTAG,ELEV,XLAPSE, 1 LNDSEA,XDATA,NSTA,KSTA,RCLOS2,LIO, 2 LS1,LS2,ERRANL(J,1),ERRANL(J,2),ERRANL(J,17), 3 ERRANL(J,18),ERRANL(J,19),IER) C IF(IER.NE.0)THEN WRITE(KFILDO,153) 153 FORMAT(' WILL NOT PROCESS THIS STATION.') C THIS FOLLOWS DIAGNOSTIC IN CLOS2. ISTOP=ISTOP+1 GO TO 130 C THIS TRANSFER MEANS THAT TWO CLOSEST STATIONS HAVE C NOT BEEN FOUND WITHIN RCLOS2 = RELVAR(1)*LIOVAR, C SO ANOTHER RANDOM STATION MUST BE USED. ENDIF C C FOR ANY STATION MORE THAN NWITH, SET FIRST DISTANCE C EQUAL TO ZERO, AND THE SECOND TO THE FIRST DISTANCE. C ALSO MODIFY OTHER ELEMENTS (COLS 17-19) TO AGREE. C THIS WON'T HAPPEN WHEN IEXTRA = 0, WHICH IS NORMAL C IF(J.GT.NWITH)THEN ERRANL(J,2)=ERRANL(J,1) ERRANL(J,1)=0. ERRANL(J,17)=0. ERRANL(J,18)=0. ERRANL(J,19)=0. ENDIF C C SET THE 2ND WORD ID FOR CLOSE STATIONS. C ID2(1)=01970000+NINT(RCLOS2) ID2(2)=02970000+NINT(RCLOS2) ID2(17)=17970000+NINT(RCLOS2) ID2(18)=18970000+NINT(RCLOS2) ID2(19)=19970000+NINT(RCLOS2) ID2(20)=20970000 C ID2(20) IS THE PREDICTAND; UUUU = 0. C D WRITE(KFILDO,154)J,CCALL(KSTA),CCALL(LS1),ERRANL(J,1), D 1 CCALL(LS2),ERRANL(J,2) D154 FORMAT(' AT 154 IN WTHOL2--', D 1 'J,CCALL(KSTA),CCALL(LS1),ERRANL(J,1),', D 2 'CCALL(LS2),ERRANL(J,2)', D 3 I3,1X,A8,4X,A8,F6.2,4X,A8,F6.2) C C COMPUTE THE VARIABILITY USING LAPSE RATE BUT NO C DISTANCE WEIGHTING. C DO 1547 M=1,2 CALL VARIW(KFILDO,CCALL,XP,YP,LTAG,LNDSEA,ELEV,XLAPSE, 1 XDATA,NSTA,KSTA,RELVAR(M)*VARTAB(LIOVAR), 2 ERRANL(J,M+2),LIO,IER) C THIS USES ONLY THE FIRST TWO VALUES IN RELVAR( ) C THERE IS CURRENTLY NO NON-ZERO ERROR RETURN. 1547 CONTINUE C C SET VALUE TO NEXT LARGER RADIUS WHEN MISSING. THERE C ARE ONLY 2 VALUES. C DO 1548 M=4,4 C IF(ERRANL(J,M).GT.9998.9)THEN ERRANL(J,M)=ERRANL(J,M-1) C GIVEN THAT THE LARGER RADII OF SEARCH ARE FOR LOWER M, C THIS BUMPS THE VALUES UP SO THAT THERE ARE NO MISSINGS. C NOTE THAT THIS IS KEYED TO SPECIFIC COLUMNS OF ERRANL, C AS THE CALL TO ELVARI. ENDIF C 1548 CONTINUE C C SET THE 2ND WORD ID FOR STATION ROUGHNESS. C ID2(3)=03970000+NINT(RELVAR(1)*VARTAB(LIOVAR)) ID2(4)=04970000+NINT(RELVAR(2)*VARTAB(LIOVAR)) C D WRITE(KFILDO,1549)(ERRANL(J,M),M=3,4) D1549 FORMAT(/,' AT 1549--(ERRANL(J,M),M=3,4)',4F10.2) C C COMPUTE A ROUGHNESS VARIABLE FOR EACH OF FOUR RADII C CENTERED ON THE GRIDPOINT CLOSEST TO THE STATION BEING C PROCESSED. THE COMPUTATION IS OVER NROUGH( ). THE C TERRAIN IN TELEV( , ) IS NOT INTERPOLATED TO THE C STATION AND EVEN GRIDPOINT INCREMENTS FROM IT. RATHER C THE GRIDPOINT VALUES IN TELEV( , ) ARE USED AS INCREMENTS C FROM THE CLOSEST GRIDPOINT TO THE STATION. HENCE, THE C NINT FUNCTION. C IXX=NINT(XP(KSTA)) JYY=NINT(YP(KSTA)) C DO 155 M=1,4 CALL ROUGH(KFILDO,CCALL(KSTA),SEALND,TELEV,NXE,NYE,IXX,JYY, 1 NROUGH(M),ERRANL(J,M+4),LIO,IER) C THERE IS NO NON-ZERO ERROR RETURN. 155 CONTINUE C C SET THE 2ND WORD ID FOR GRIDPOINT ROUGHNESS. C ID2(5)=05970000+NROUGH(1) ID2(6)=06970000+NROUGH(2) ID2(7)=07970000+NROUGH(3) ID2(8)=08970000+NROUGH(4) C C COMPUTE THE FOUR DATA VARIABILITY VARIABLES WITH C NO LAPSE OR DISTANCE WEIGHTING. THE SEARCH IS OVER RELVAR( ). C DO 1552 M=1,4 CALL VARI(KFILDO,CCALL,XP,YP,LTAG,LNDSEA,XDATA,NSTA,KSTA, 1 RELVAR(M)*VARTAB(LIOVAR),ERRANL(J,M+8),LIO,IER) C THERE IS CURRENTLY NO NON-ZERO ERROR RETURN. 1552 CONTINUE C C SET VALUE TO NEXT LARGER RADIUS WHEN MISSING. C DO 1554 M=10,12 C IF(ERRANL(J,M).GT.9998.9)THEN ERRANL(J,M)=ERRANL(J,M-1) C GIVEN THAT THE LARGER RADII OF SEARCH ARE FOR LOWER M, C THIS BUMPS THE VALUES UP SO THAT THERE ARE NO MISSINGS. C NOTE THAT THIS IS KEYED TO SPECIFIC COLUMNS OF ERRANL, C AS THE CALL TO VARI. ENDIF C 1554 CONTINUE C C SET THE 2ND WORD ID FOR DATA VARIABILITY. C ID2(09)=09970000+NINT(RELVAR(1)*VARTAB(LIOVAR)) ID2(10)=10970000+NINT(RELVAR(2)*VARTAB(LIOVAR)) ID2(11)=11970000+NINT(RELVAR(3)*VARTAB(LIOVAR)) ID2(12)=12970000+NINT(RELVAR(4)*VARTAB(LIOVAR)) C D WRITE(KFILDO,1555)(ERRANL(J,M),M=9,12) D1555 FORMAT(/,' AT 1555--(ERRANL(J,M),M=9,12)',4F10.2) C C COMPUTE THE VARIABILITY USING WEIGHTING FUNCTION AND LAPSE. C DO 1557 M=1,4 CALL VARIL(KFILDO,CCALL,XP,YP,LTAG,LNDSEA,ELEV,XLAPSE, 1 XDATA,NSTA,KSTA,RELVAR(M)*VARTAB(LIOVAR), 2 ERRANL(J,M+12),LIO,IER) C THERE IS CURRENTLY NO NON-ZERO ERROR RETURN. 1557 CONTINUE C C SET VALUE TO NEXT LARGER RADIUS WHEN MISSING. C DO 1558 M=14,16 C IF(ERRANL(J,M).GT.9998.9)THEN ERRANL(J,M)=ERRANL(J,M-1) C GIVEN THAT THE LARGER RADII OF SEARCH ARE FOR LOWER M, C THIS BUMPS THE VALUES UP SO THAT THERE ARE NO MISSINGS. C NOTE THAT THIS IS KEYED TO SPECIFIC COLUMNS OF ERRANL, C AS THE CALL TO ELVARI. ENDIF C 1558 CONTINUE C C SET THE 2ND WORD ID FOR STATION ROUGHNESS. C ID2(13)=13970000+NINT(RELVAR(1)*VARTAB(LIOVAR)) ID2(14)=14970000+NINT(RELVAR(2)*VARTAB(LIOVAR)) ID2(15)=15970000+NINT(RELVAR(3)*VARTAB(LIOVAR)) ID2(16)=16970000+NINT(RELVAR(4)*VARTAB(LIOVAR)) C D WRITE(KFILDO,1559)(ERRANL(J,M),M=13,16) D1559 FORMAT(/,' AT 1559--(ERRANL(J,M),M=13,16)',4F10.2) C WSTA(J)=CCALL(KSTA) WRITE(DUMCAL(J)(1:8),156,IOSTAT=IOS,ERR=157)J 156 FORMAT('DUM',I2.2,' ') CCC WRITE(KFILDO,1565)J,KSTA,DUMCAL(J) CCC 1565 FORMAT(' AT 1656--J,KSTA,DUMCAL(J)',I4,I6,2X,A8) GO TO 160 C 157 WRITE(KFILDO,158)IOS 158 FORMAT(/,' ****ERROR ON INTERNAL WRITE IN WTHOL2 AT 158', 1 ' IOSTAT=',I5,'. PROCEEDING.') ISTOP=ISTOP+1 C C AT THIS POINT, FOR J = NWITH + IEXTRA: C ERRANL(J,1) AND ERRANL(J,2) HOLD THE DISTANCES TO THE C TWO CLOSEST STATIONS. C ERRANL(J,3) AND ERRANL(J,4) ARE NOT USED. C ERRANL(J,5) THROUGH ERRANL(J,8) HOLD THE FOUR C ROUGHNESS VARIABLES. C ERRANL(J,9) THROUGH ERRANL(J,12) HOLD THE DATA C VARIABILITY VARIABLES. C ERRANL(J,13) THROUGH ERRANL(J,16) HOLD THE ELEVATION C VARIABILITY VARIABLES. C ERRANL(J,17) IS THE ABSOLUTE DIFFERENCE IN ELEVATION C BETWEEN THE STATION AND ITS CLOSEST NEIGHBOR. C ERRANL(J,18) IS THE ABSOLUTE DIFFERENCE BETWEEN THE C STATION DATA VALUE AND VALUE OF ITS NEIGHBOR AFTER C THE LAPSE HAS BEEN APPLIED. C ERRANL(J,19) IS THE PRODUCT OF ERRANL(J,17) AND C ERRANL(J,1). C ERRANL(J,20) WILL BE THE PREDICTAND. C THE FIRST NWITH ARE WITHHELD. NOTE THAT THEY HAVE NOT C BEEN WITHHELD FOR THE ABOVE CALLS. C 160 WHOLD(KSTA)=XDATA(KSTA) C ALL NWITH+IEXTRA STATIONS ARE PUT INTO WHOLD( ). THIS C IS FOR EASE OF CALCULATIONS IN SUBROUTINE DIFWH. C IF(J.LE.NWITH)THEN C THIS IS A POINT TO WITHHOLD. NOTE THE IEXTRA POINTS C ARE PUT INTO THE LIST ABOVE BUT ARE NOT WITHHELD. C XDATA(KSTA)=9999. LTAGWH(KSTA)=LTAG(KSTA) C NOTE THAT LTAGWH( ) IS NOT INITIALIZED. SHOULD NOT HAVE C TO BE. LTAG(KSTA)=1 D WRITE(KFILDO,168)KSTA,XDATA(KSTA),WHOLD(KSTA),LTAGWH(KSTA), D 1 LTAG(KSTA),XP(KSTA),YP(KSTA) D168 FORMAT(/,' AT 168 IN WTHOL2--KSTA,XDATA(KSTA),WHOLD(KSTA),', D 1 'LTAGWH(KSTA),LTAG(KSTA),XP(KSTA),YP(KSTA)',/, D 2 I6,2F10.2,2I7,2F9.2) D ELSE D WRITE(KFILDO,1680)J,KSTA,XDATA(KSTA),WHOLD(KSTA), D 1 XP(KSTA),YP(KSTA) D1680 FORMAT(/,' AT 1680 IN WTHOL2--J,KSTA,XDATA(KSTA),', D 1 'WHOLD(KSTA),XP(KSTA),YP(KSTA)',/,2I6,4F10.2) ENDIF C D WRITE(KFILDO,169)CCALL(KSTA) D169 FORMAT(/,' WITHHELD OR EXTRA STATION = ',A8) C ENDIF C 180 CONTINUE C IF(MFIRST.NE.0)GO TO 200 C C THIS SECTION DELETES withheld IF IT EXISTS, REOPENS IT, C AND WRITES THE CALL LETTERS RECORD. IP24 IS USED SO C THAT A NEW OUTPUT FILE WOULD NOT BE NECESSARY IN THE C U155 CONTROL FILE FOR THIS SPECIAL STUDY. C FLNAM(1:8)='withheld' CALL SYSTEM('rm -f '//FLNAM) C REMOVE AN EXISTING FILE OF NAME FLNAM IF IT EXISTS. C AN ERROR WILL OCCUR IF THE FILE ALREADY EXISTS. D WRITE(KFILDO,185)IP24 D185 FORMAT(/,' AT 185 IN WTHOL2--IP24,FLNAM',I10,2X,A8) COPS OPEN(UNIT=IP24,FILE=FLNAM,STATUS='NEW',FORM='UNFORMATTED', OPEN(UNIT=IP24,STATUS='NEW',FORM='UNFORMATTED', 1 IOSTAT=IOS,ERR=900) C C USE WRDIR TO WRITE THE CALL LETTERS RECORD ON FTN24. C KWRITE=1 NTOTBY=0 NTOTRC=0 CALL WRDIR(KFILDO,IP24,KWRITE,DUMCAL,NWITH+IEXTRA, 1 NTOTBY,NTOTRC,L3264W,IER) C IF(IER.NE.0)THEN WRITE(KFILDO,190) 190 FORMAT(/,' ****FILE WITHHELD COULD NOT BE INITIALIZED', 1 ' WITH CALL LETTERS. FATAL ERROR IN WTHOL2.') CALL W3TAGE('WTHOL2') STOP 190 ENDIF C MFIRST=1 C 200 CONTINUE C STATEMENT 200 PROVIDES FOR THE EXECUTION OF TIMPR C ON OTHER THAN THE FIRST ENTRY (WHEN COMPILED WITH C /D). C D WRITE(KFILDO,298)(WSTA(J),(ERRANL(J,M),M=1,20), D 1 J=1,NWITH+IEXTRA) D298 FORMAT(/' STATIONS SELECTED AND WITHHELD DATA',/, D 1 (A8,18F6.1,2F8.1)) WRITE(KFILDO,299)NDATE,(WSTA(J),J=1,NWITH+IEXTRA) 299 FORMAT(/,' STATIONS SELECTED IN WHTOL2 FOR DATE/TIME',I12,/, 1 (10(2X,A8))) CALL TIMPR(KFILDO,KFILDO,'END WTHOL2 ') RETURN C 900 WRITE(KFILDO,901)FLNAM 901 FORMAT(/,' ****ERROR OPENING FILE ',A8,' ON UNIT FTN24', 1 ' IN WTHOL2. FATAL ERROR.') CALL W3TAGE('WTHOL2') STOP 900 END