SUBROUTINE WTHOL3(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 SEPTEMBER 2009 IM INSERTED CALL TO SYSTEM CLOCK
C MARCH 2010 SDS MODIFIED OPEN STATEMENT FOR OPERATIONS
C JUNE 2014 GLAHN INSERTED IER =0
C JUNE 2014 GLAHN INSERTED CALL TO W3TAGE BEFORE STOPS
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 THIS WHTOL3 IS THE SAME AS WTHOL2 EXCEPT THE FIRST SEED
C TO USE COMES FROM THE SYSTEM CLOCK RATHER THAN FROM THE
C INPUT VARIABLE SEED. IF WHTOL3 IS ENTERED ONLY ONCE
C FOR A RUN, THE ONLY SEED USED COMES FROM THE SYSTEM CLOCK.
C IF IT IS USED MORE THAN ONCE, THE SEED IS SAVED AS IN
C WTHOL2.
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 (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 NEIGHBOR AFTER LAPSE
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
C IER=0
C
WRITE(KFILDO,105)NROUGH,RELVAR,VARTAB
105 FORMAT(' AT 105 IN WTHOL3--NROUGH,RELVAR,VARTAB',/,
1 4I4,6F10.1)
C
CALL TIMPR(KFILDO,KFILDO,'START WTHOL3 ')
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)
IF(JFIRST.EQ.0)THEN
C
C USE SYSTEM_CLOCK TO GET RANDOM SEED.
C MAXIMUM VALUE FOR ICOUNT IS 8639999 FOR CURRENT IBM SYSTEM
C (I.E., CIRRUS/STRATUS), SO ICOUNT WILL FIT FOR SEED
C (4 BYTE REAL VALUE)
C
CALL SYSTEM_CLOCK(ICOUNT)
SEED=FLOAT(ICOUNT)
CALL SRAND(SEED)
C THIS SETS THE SEED FOR THE RANDOM NUMBER GENERATOR.
C
C BECAUSE THE SYSTEM CLOCK SEED MAY HAVE A PATTERN
C (ONE LOOKS A LOT LIKE THE PREVIOUS ONE), THE LOOP BELOW
C IS TO MAKE SURE THE SEED USED IS "MORE RANDOM." THIS LOOP
C IS ALSO USED IN WTHOL2, AND HAS THE SAME PURPOSE.
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 WTHOL3--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 WTHOL3--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 WTHOL3.')
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 WTHOL3',
2 'AT 146.')
CALL W3TAGE('WTHOL3')
STOP 146
C
ELSE
D WRITE(KFILDO,147)X,Y,CPNDFD(NINT(X),NINT(Y))
D147 FORMAT(' IN WTHOL3 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 WTHOL3',
2 'AT 1475.')
CALL W3TAGE('WTHOL3')
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 WTHOL3--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 WTHOL3--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
!STOP 11111
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 WTHOL3--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 WTHOL3--',
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 WTHOL3 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 WTHOL3--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 WTHOL3--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 WTHOL3--IP24,FLNAM',I10,2X,A8)
COPS OPEN(UNIT=IP24,FILE=FLNAM,STATUS='NEW',FORM='UNFORMATTED',
COPS 1 IOSTAT=IOS,ERR=900)
OPEN(UNIT=IP24,STATUS='NEW',FORM='UNFORMATTED',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 WTHOL3.')
CALL W3TAGE('WTHOL3')
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 WTHOL3 ')
RETURN
C
900 WRITE(KFILDO,901)FLNAM
901 FORMAT(/,' ****ERROR OPENING FILE ',A8,' ON UNIT FTN24',
1 ' IN WTHOL3. FATAL ERROR.')
CALL W3TAGE('WTHOL3')
STOP 900
END