SUBROUTINE POINTS(KFILDO,CCALL,NAME,ELEV,IWBAN,STALAT,STALON, 1 LNDSEA,IQUAL,NOPAR,ND1,NSTA,ISMPL, 2 TELEV,SEALND,NXE,NYE,BMESH,EMESH,NXL,NYL, 3 NPROJ,ORIENT,XLAT,ALATL,ALONL, 4 ISTOP,IER) C C OCTOBER 2005 GLAHN TDL MOS-2000 C OCTOBER 2005 GLAHN ADDED NOPAR( ) C MAY 2007 GLAHN INCREASED IQUAL( ,2) TO IQUAL( ,5) C DECEMBER 2008 GLAHN HP FUNCTION RAN REPLACED WITH IBM C 'RANDOM_NUMBER' AND 'RANDOM_SEED' C FEBRUARY 2009 GLAHN REPLACED THE IBM RANDOM NUMBER C GENERATOR WITH ANOTHER SIMPLER ONE C *****HASN'T BEEN TESTED**** C C PURPOSE C TO GENERATE DATA POINTS TO ANALYZE FROM A FIRST GUESS GRID. C THIS CAN ADD TO A LIST OF STATION DATA ALREADY READ C IN IF NSTA GT 0. ARRAYS CCALL( , ), ELEV( ), IWBAN( ), C STALAT( ), STALON( ), LNDSEA( ), AND IQUAL( ) ARE C FILLED. THIS ROUTINE CALLED FROM U155. BECAUSE SEEDT C IS SET TO THE SAME VALUE FOR EACH RUN, THE POINTS C GENERATED WILL BE THE SAME FOR EACH ANALYSIS FOR WHICH C THEY ARE USED. THIS IS PROBABLY AN ADVANTAGE FOR C CONSISTENCY. 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 CCALL(K,J) = 8-CHARACTER STATION CALL LETTERS (OR GRIDPOINT C LOCATIONS FOR GRID DEVELOPMENT) TO PROVIDE C OUTPUT FOR (J=1) AND 5 POSSIBLE OTHER STATION C CALL LETTERS (J=2,6) THAT CAN BE USED INSTEAD C IF THE PRIMARY (J=1) STATION CANNOT BE FOUND C IN AN INPUT DIRECTORY (K=1,NSTA). ALL STATION C DATA ARE KEYED TO THIS LIST. (CHARACTER*8) C (INPUT/OUTPUT) C NAME(K) = NAMES OF STATIONS (K=1,NSTA). (CHARACTER*20) C (INPUT/OUTPUT) C ELEV(K) = ELEVATION OF STATIONS IN METERS (K=1,NSTA). C (INPUT/OUTPUT) C IWBAN(K) = WBAN NUMBERS OF STATIONS (K=1,NSTA). C (INPUT/OUTPUT) C STALAT(K) = LATITUDE OF STATIONS (K=1,NSTA). (INPUT/OUTPUT) C STALON(K) = LONGITUDE OF STATIONS (K=1,NSTA). 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) GRIDPONTS. C (INPUT/OUTPUT) C IQUAL(K,I) = THE QUALITY VALUES FROM THE STATION DICTIONARY C FOR FIVE POSSIBLE DATA TYPES (K=1,ND1) (I=1,5). C SET = 4. (OUTPUT) C NOPAR(K) = THE NUMBER OF PAIRS FOR STATION K. FOR SAMPLED C POINTS, THIS IS SET TO 9999. (INPUT/OUTPUT) C ND1 = MAXIMUM NUMBER OF STATIONS THAT CAN BE DEALT C WITH. (INPUT) C NSTA = THE NUMBER OF STATIONS BEING DEALT WITH. THE C NUMBER OF VALUES IN CCALL( , ), ETC. ON INPUT, C THIS IS THE NUMBER OF STATIONS IN THE DIRECTORY C READ; ON OUTPUT, IT IS NSTA+ISMPL. C (INPUT/OUTPUT) C ISMPL = THE NUMBER OF DATA POINTS TO ADD. (INPUT) C TELEV(IX,JY) = THE TERRAIN ELEVATION FROM THE MOS-2000 EXTERNAL C RANDOM ACCESS FILE AT MESHLENGTH EMESH C (IX=1,NXE) (JY=1,NYE). (INPUT) C SEALND(IX,JY) = THE LAND/SEA MASK AT MESHLENGTH EMESH C (IX=1,NXE) (JY=1,NYE)). C 0 = OCEAN WATER GRIDPOINTS; C 3 = INLAND WATER GRIDPOINTS. C 9 = LAND GRIDPOINTS. C (INPUT) C NXE = X-EXTENT OF TELEV( , )) AND SEALND( , ). C (INPUT) C NYE = Y-EXTENT OF TELEV( , )) AND SEALND( , ). C (INPUT) C BMESH = ACTUAL MESH LENGTH IN KM CORRESPONDING TO C MESHB. (INPUT) C EMESH = ACTUAL MESH LENGTH IN KM CORRESPONDING TO C MESHE. (INPUT) C NPROJ = MAP PROJECTION. (INPUT) C ORIENT = ORIENTATION OF GRID IN WEST LONGITUDE. (INPUT) C XLAT = NORTH LATITUDE AT WHICH GRIDLENGTH IS SPECIFIED. C (INPUT) C ALATL = NORTH LATITUDE OF LOWER LEFT CORNER POINT C OF A GRID OF THE SIZE NXL, NYL. TRUNCATED C TO TEN THOUSANDTHS OF DEGREES. NOTE THAT THE C MOS-2000 ARCHIVE IS ONLY TO THOUSANDTHS OF C DEGREES. (INPUT) C ALONL = WEST LONGITUDE OF LOWER LEFT CORNER POINT C OF A GRID OF THE SIZE NXL, NYL. TRUNCATED C TO TEN THOUSANDTHS OF DEGREES. NOTE THAT THE C MOS-2000 ARCHIVE IS ONLY TO THOUSANDTHS OF C DEGREES. (INPUT) C NXL = THE SIZE OF THE ANALYSIS GRID FOR THIS RUN C IN THE X DIRECTION IN MESHB UNITS. (INPUT) C NYL = THE SIZE OF THE ANALYSIS GRID FOR THIS RUN C IN THE Y DIRECTION IN MESHB UNITS. (INPUT) C ISTOP = 0 MEANS THE PROGRAM IS RUNNING OK UP TO THIS C POINT. WHENEVER AN ERROR OCCURS THAT SHOULD C HALT THE PROGRAM AFTER INPUT DIAGNOSTICS ARE C PRINTED, ISTOP IS SET = ISTOP+1. (INPUT/OUTPUT) C IER = STATUS RETURN. C 0 = GOOD RETURN. C 777 = FATAL ERROR. C OTHER VALUES CAN COME FROM CALLED SUBROUTINES. C (OUTPUT) C RMESH = RATIO OF MESH TO MESHE. (INTERNAL) C SEED = A REAL NUMBER THAT STARTS THE RANDOM NUMBER C SEQUENCE. (INTERNAL) C IWLOOP = THE NUMBER OF TIMES TO CALL THE RANDOM NUMBER C GENERATOR BEFORE THE VALUES ARE USED. THIS IS C TO GET BEYOND ANY POSSIBLE FLAW IN THE VALUE C OF ISEED. SET BY DATA STATEMENT. (INTERNAL) C C NONSYSTEM SUBROUTINES USED C LMIJLL, PSIJLL, MCIJLL C 1 2 3 4 5 6 7 X C CHARACTER*8 CCALL(ND1,6) CHARACTER*20 NAME(ND1) C DIMENSION ELEV(ND1),IWBAN(ND1),STALAT(ND1),STALON(ND1), 1 IQUAL(ND1,5),LNDSEA(ND1),NOPAR(ND1) DIMENSION TELEV(NXE,NYE),SEALND(NXE,NYE) C DATA SEED/.9192225/ C C X IS MODIFIED AND SAVED FOR THE NEXT ITERATION. DATA ICOUNT/0/ C ICOUNT COUNTS THE NUMBER OF ITEMS ADDED. DATA IWLOOP/10/ C IWLOOP IS THE NUMBER OF TIMES TO CALL THE RANDOM NUMBER C GENERATOR BEFORE THE VALUES ARE USED. CALL SRAND(SEED) C THIS SETS THE SEED FOR THE RANDOM NUMBER GENERATOR. C RMESH=REAL(NINT(BMESH/EMESH)) C RMESH IS THE RATIO OF THE MESH LENGTH OF THE ANALYSIS GRID C C TO THE TERRAIN GRID. C CALL RAND IWLOOP TIMES TO GET PAST ANY PROBLEMS WITH USING AN C INAPPROPRIATE INITIAL ISEED. C DO 120 J=1,IWLOOP X=RAND() 120 CONTINUE C CD WRITE(KFILDO,121)ISEED,X,RMESH,NSTA,ISMPL,NXL,NYL CD121 FORMAT(/' AT 121 IN POINTS--ISEED,X,RMESH,NSTA,ISMPL,NXL,NYL', CD 1 I12,2F10.6,4I6) C DATA POINTS ARE GENERATED BY A RANDOM SAMPLING FROM THE GRID C OVER WHICH THE ANALYSIS IS BEING DONE. C DO 200 K=NSTA+1,NSTA+ISMPL C IF(K.GT.ND1)THEN WRITE(KFILDO,125)ND1,NSTA,ISMPL 125 FORMAT(/' *****ND1 ABOUT TO BE EXCEEDED IN POINTS. ND1 =',I6, 1 ' ISMPL =',I6,'. PROCEEDING.') NSTA=K-1 ISTOP=ISTOP+1 GO TO 600 ENDIF C NOPAR(K)=9999 C NOPAR( ) SIGNIFIES THERE IS NO PAIR LIST FOR THIS POINT. ICOUNT=ICOUNT+1 130 X=RAND() C X VARIES FROM 0 TO 1. X=X*NXL C X IS THE LOCATION IN THE GRID IN THE X DIRECTION. CD WRITE(KFILDO,132)K,ISEED,X CD132 FORMAT(' AT 132 IN POINTS--K,ISEED,X',I6,I12,F12.4) IF(X.LT.1..OR.X.GE.NXL)GO TO 130 C THIS IS FOR SAFETY. ZERO COULD OCCUR. C C PICK THE Y POSITION. C 140 Y=RAND() Y=Y*NYL C Y IS THE LOCATION IN THE GRID IN THE Y DIRECTION. CD WRITE(KFILDO,142)K,ISEED,Y CD142 FORMAT(' AT 142 IN POINTS--K,ISEED,Y',I6,I12,F12.4) IF(Y.LT.1..OR.Y.GE.NYL)GO TO 140 C THIS IS FOR SAFETY. ZERO COULD OCCUR. C C FILL CALL LETTERS ARRAY CCALL( , ). USE NUMBERS STARTING C WITH 'Q_'. C WRITE(CCALL(K,1)(1:8),145,IOSTAT=IOS,ERR=900)ICOUNT 145 FORMAT('Q_',I6.6) CCALL(K,2)=' ' CCALL(K,3)=' ' CCALL(K,4)=' ' CCALL(K,5)=' ' CCALL(K,6)=' ' C C FILL THE TWO QUALITY COLUMNS IQUAL( , ) AND IWBAN( ). C IQUAL(K,1)=4 IQUAL(K,2)=4 IQUAL(K,3)=4 IQUAL(K,4)=4 IQUAL(K,5)=4 IWBAN(K)=0 C C CALCULATE THE LATITUDE AND LONGITUDE OF THE POINT. C IF(NPROJ.EQ.3)THEN CALL LMIJLL(KFILDO,X,Y, 1 BMESH*1000.,ORIENT,XLAT, 2 ALATL,ALONL,ALATP,ALONP,IER) CD WRITE(KFILDO,147)BMESH,ORIENT,XLAT,ALATL,ALONL,ALATP,ALONP CD147 FORMAT(' IN POINTS AT 147--BMESH,ORIENT,XLAT,ALATL,', CD 1 'ALONL,ALATP,ALONP',7F10.2) C IF(IER.NE.0)THEN WRITE(KFILDO,150)IER 150 FORMAT(/' FATAL ERROR IN LMIJLL FROM POINTS,', 1 ' AT 150, IER =',I4) ISTOP=ISTOP+1 IER=777 GO TO 600 ENDIF C ELSEIF(NPROJ.EQ.5)THEN CALL PSIJLL(KFILDO,X,Y, 1 BMESH*1000.,ORIENT,XLAT, 2 ALATL,ALONL,ALATP,ALONP) ELSEIF(NPROJ.EQ.7)THEN CALL MCIJLL(KFILDO,X,Y, 1 BMESH*1000.,XLAT, 2 ALATL,ALONL,ALATP,ALONP) ELSE WRITE(KFILDO,151)NPROJ 151 FORMAT(/' ****MAP PROJECTION NUMBER NPROJ NOT', 1 ' 3, 5, OR 7. FATAL ERROR IN POINTS AT 151.') ISTOP=ISTOP+1 IER=777 GO TO 600 ENDIF C C FILL THE STALAT( ) AND STALON( ) ARRAYS. C STALAT(K)=ALATP STALON(K)=ALONP C C FILL NAME ARRAY NAME( ). USE LATITUDE AND LONGITUDE *10. C ILAT=NINT(ALATP*10.) ILON=NINT(ALONP*10.) WRITE(NAME(K)(1:20),155,IOSTAT=IOS,ERR=900)ILAT,ILON 155 FORMAT('LAT ',I4.4,'; LON ',I4.4,' ') C C PUT THE ELEVATION OF THE POINT INTO ELEV( ). C INTERPOLATE LINEARLY. C IF(RMESH.EQ.1.)THEN XE=X YE=Y ELSE XE=(X-1.)*RMESH+1. YE=(Y-1.)*RMESH+1. ENDIF CALL ITRP(TELEV,NXE,NYE,XE,YE,ELEV(K)) CD WRITE(KFILDO,157)K,NXE,NYE,X,Y,ELEV(K) CD157 FORMAT(' AT 156 IN POINTS--K,NXE,NYE,X,Y,ELEV(K)', CD 1 3I6,3F10.3) C C PUT THE LAND/SEA TYPE INTO LNDSEA( ). C TAKE THE CLOSEST POINT. C IX=NINT(XE) JY=NINT(YE) C LNDSEA(K)=SEALND(IX,JY) CD WRITE(KFILDO,158)K,NXE,NYE,X,Y,LNDSEA(K) CD158 FORMAT(' AT 156 IN POINTS--K,NXE,NYE,X,Y,LNDSEA(K)', CD 1 3I6,2F10.3,I4) 200 CONTINUE C NSTA=NSTA+ISMPL 600 RETURN C C 900 WRITE(KFILDO,901)IOS 901 FORMAT(/,' ****ERROR ON INTERNAL WRITE IN POINTS AT 900', 1 ' IOSTAT=',I5,'. FATAL ERROR.') ISTOP=ISTOP+1 IER=777 RETURN END