SUBROUTINE WTHOL1(KFILDO,IP24,ID,NDATE,CCALL,XDATA,ND1, 1 LTAG,LNDSEA,WHOLD,LTAGWH,XP,YP,NSTA,NX,NY,MESH, 2 CPNDFD,NXE,NYE,MESHE,NCLIPY, 3 SEED,NWITH,ISTOP,IER) C C JUNE 2005 GLAHN TDL MOS-2000 C JULY 2005 GLAHN ADDED CCALL TO CALL AND KFIL8 C AUGUST 2005 GLAHN ADDED IWLOOP TO CALL C AUGUST 2005 GLAHN CHANGED LNDSEA(K).NE.4 TO C LNDSEA(K).LE.3 C APRIL 2006 GLAHN ADDED IP24 TO CALL; ELIMINATED KFIL08 C OCTOBER 2006 GLAHN CHANGED LNDSEA(K).LE.3 TO LE.6; C COMMENT FOR LNDSEA( ) C OCTOBER 2006 GLAHN ADDED MESH,CPNDFD,NXE,NYE,MESHE, C NCLIPO,ID,NDATE; COMMENTS C NOVEMBER 2006 GLAHN MODIFIED TO ONLY WITHHOLD STATIONS C WITH IDENTIFIERS STARTING WITH ZK C DECEMBER 2006 GLAHN ADDED NOTE UNDER PURPOSE C JANUARY 2007 GLAHN REVISED NOTE UNDER PURPOSE C JANUARY 2007 GLAHN ADDED ND1 TO CALL C FEBRUARY 2007 GLAHN ADDED CHECK FOR LTAG(K) NE 0 IN C DO 150 LOOP C FEBRUARY 2007 GLAHN REMOVED DEPENDENCY ON "Z" C JUNE 2007 GLAHN MOVED LTAG TO LINE 2 IN CALL C SEPTEMBER 2007 GLAHN FORMAT CHANGE AT 146 C NOVEMBER 2007 GLAHN REMOVED IWLOOP; CHANGED NAME FROM C WITHOL TO WTHOL1 C DECEMBER 2007 GLAHN SPELL CHECK 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' AND 'RANDOM_SEED' C JANUARY 2009 GLAHN RANDOM NUMBER GENERATOR NOW RAND( ) C FEBRUARY 2009 GLAHN CALLED SRAND(SEED) ONLY WHEN JFIRST =0 C JUNE 2014 GLAHN IINSERTED IER=0 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: THE CURRENT DATA ARE ALWAYS PUT INTO XDATA( ,2). C THE DATA TO ANALYZE ARE MERGED INTO XDATA( ,1). C NOTE: FOR VERIFYING ACCURACY OF FORECASTS AND THE MOST C ACCURATE REPORTING STATIONS ARE DESIRED, DO NOT C WITHHOLD STATIONS UNLESS CCALL( ) STARTS WITH "K". 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. (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 XDATA(K,J) = ON INPUT, HOLDS THE FULL SET OF DATA (K=1,NSTA) C (J=1,2). THE FIRST COLUMN IS THE DATA TO C ANALYZE; THE SECOND COLUMN IS THE ON-CYCLE C DATA. IF THERE IS ONLY ONE CYCLE INVOLVED, C COLUMNS 1 AND 2 ARE THE SAME. C ON OUTPUT, THE WITHHELD DATA HAVE BEEN SET TO C MISSING = 9999 IN XDATA( ,1). (INPUT/OUTPUT) 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 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 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 NSTA = NUMBER OF STATIONS OR LOCATIONS BEING DEALT C WITH. (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 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 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 SEED = SEED FOR THE RANDOM NUMBER GENERATOR. C (INPUT/OUTPUT) C NWITH = NUMBER OF STATIONS TO WITHHOLD. (INPUT) C IER = 0 = GOOD RETURN. (OUTPUT) C ISTOP = INCREMENTED BY 1 WHEN AN ERROR OCCURS. C (INPUT/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. C JFIRST = 0 ON FIRST ENTRY. AFTER RAN IS CALLED 10 TIMES C IT IS SET TO 1. (INTERNAL) (SAVED) C 1 2 3 4 5 6 7 X C C NONSYSTEM SUBROUTINES CALLED C NONE. C CHARACTER*1 Q CHARACTER*8 CCALL(NSTA) C DIMENSION ID(4) DIMENSION XDATA(ND1,2),LTAG(NSTA),LNDSEA(NSTA),WHOLD(NSTA), 1 LTAGWH(NSTA),XP(NSTA),YP(NSTA) DIMENSION CPNDFD(NXE,NYE) C DATA IFIRST/0/, 1 JFIRST/0/ C SAVE IFIRST,JFIRST C D CALL TIMPR(KFILDO,KFILDO,'START WTHOL1 ') C IER=0 C IF(IP24.NE.0)THEN WRITE(IP24,110)NWITH,(ID(J),J=1,4),NDATE 110 FORMAT(/,I4,' STATIONS ARE BEING WITHHELD PER ANALYSIS FOR', 1 ' VARIABLE',3I10.9,I10.3,' FOR DATE',I12) ENDIF 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 DO 180 J=1,NWITH 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 WITHOL--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 WITHOL--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 NOTE THAT THIS GENERATES ANOTHER X VALUE IN ADDITION C TO A Y VALUE. oTHERWISE BIAS CAN OCCUR. C C THROW OUT A GRIDPOINT THAT IS NOT WITHIN THE NDFD GRID. 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 WTHOL1.') IFIRST=1 ISTOP=ISTOP+1 ENDIF C ELSEIF(MESH.NE.MESHE)THEN C IF(IFIRST.EQ.0)THEN 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 ' THIS WILL LIKELY BIAS RESULTS.', 4 ' PROCEEDING AT 146 IN WTHOL1.') IFIRST=1 ISTOP=ISTOP+1 ENDIF C ELSE D WRITE(KFILDO,147)X,Y,CPNDFD(NINT(X),NINT(Y)) D147 FORMAT(' IN WITHOL 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. ENDIF C C FIND THE CLOSEST NON-WATER STATION WITH NON-MISSING DATA C FOR THE ON-CYCLE RUN. THIS WILL GUARANTEE A MATCH WHEN C COMPARING SINGLE ON-CYCLE DATA WITH A COMBINATION OF C RUNS. C DISTSQ=99999999. KSTA=9999 C DO 150 K=1,NSTA IF(XDATA(K,2).EQ.9999.)GO TO 150 C DON'T DEAL WITH A MISSING DATUM FOR THE ON-CYCLE RUN. IF(LTAG(K).NE.0)GO TO 150 IF(LNDSEA(K).LE.6)GO TO 150 C DON'T ELIMINATE ANY POINTS THAT CAN BE USED FOR WATER. IF(XP(K).LT.1..OR.XP(K).GE.NX)GO TO 150 IF(YP(K).LT.1..OR.YP(K).GE.NY)GO TO 150 C DON'T WITHHOLD STATIONS OUTSIDE GRID. THIS INSURES C 4 GRIDPOINTS ARE AVAILABLE FOR LINEAR INTERPOLATION. C HOWEVER, IF THE POINT IS OUTSIDE THE FIRST GUESS AREA, C AN ANALYSIS WON'T BE DONE, AND THE WITHHELD STATION C CAN'T BE "VERIFIED." Q=CCALL(K) C D WRITE(KFILDO,1499)K,CCALL(K),Q D1499 FORMAT(' AT 1499 IN WTHOL1--K,CCALL(K),Q',I6,2X,A8,2X,A1) C IF(Q.NE.'K')GO TO 150 C WHEN VERIFYING ACCURACY OF MOS FORECASTS, USE ONLY THE C BEST OBS; DON'T WITHHOLD STATIONS UNLESS THE IDENTIFIER C STARTS WITH 'K'. THIS TEST WILL THROW OUT MANY OF THE C STATIONS. THE LETTER Q IS USED BECAUSE K IS OTHERWISE C USED. C WHEN VERIFYING THE FIT OF THE ANALYSIS TO THE DATA, C DON'T MAKE THIS RESTRICTION. 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,155)X,Y 155 FORMAT(/' COULD NOT FIND CLOSEST STATION TO GRIDPOINT', 1 ' LOCATION =',2F7.2,'. TRY ANOTHER GRIDPOINT.') GO TO 130 ELSE C C THIS IS A POINT TO WITHHOLD. C WHOLD(KSTA)=XDATA(KSTA,1) XDATA(KSTA,1)=9999. LTAGWH(KSTA)=LTAG(KSTA) C NOTE THAT LTAGWH( ) IS NOT INITIALIZED. SHOULD NOT HAVE C TO BE. LTAG(KSTA)=1 C D WRITE(KFILDO,158)KSTA,XDATA(KSTA,1),WHOLD(KSTA),LTAGWH(KSTA), D 1 LTAG(KSTA),XP(KSTA),YP(KSTA) D158 FORMAT(/' AT 158 IN WITHOL--KSTA,XDATA(KSTA),WHOLD(KSTA),', D 1 'LTAGWH(KSTA),LTAG(KSTA),XP(KSTA),YP(KSTA)', D 2 I6,2F10.2,2I4,2F9.2) C IF(IP24.NE.0)THEN WRITE(IP24,159)CCALL(KSTA) 159 FORMAT(A8) ENDIF C ENDIF C 180 CONTINUE C D CALL TIMPR(KFILDO,KFILDO,'END WTHOL1 ') RETURN END