SUBROUTINE SKYMBO(KFILDO,KFIL10,NDATE,ID,IDPARS,JD, 1 IDCIG,XDATA,XP,YP,FD2,NVAL, 2 NCAT,CONST,NSCAL,RAD, 3 LSTORE,ND9,LITEMS, 4 IS0,IS1,IS2,IS4,ND7, 5 IPACK,IWORK,DATA,ND5, 6 CORE,ND10,NBLOCK,NFETCH, 7 L3264B,ISTOP,IER) C C APRIL 2011 GLAHN MDL MOS-2000 C JUNE 2011 IM MODIFIED COMMENTS FROM C CATEGORICAL CEILING HEIGHT C TO CEILING HEIGHT IN HUNDREDS C OF FEET C OCTOBER 2011 G. WAGNER ADJUSTED CALL TO GFETCH FOR C CEILING HEIGHTS TO USE IERCIG FOR C ERROR HANDLING, TO PREVENT A C FATAL IER VALUE UPON RETURN TO C U405A. AT REQUEST OF LAMP, C ADDED A WARNING WHEN NO CIG ID C IS SPECIFIED IN THE U405A CN. C FEBRUARY 2013 IM ADDED '.AND.XDATA(L).LT.9998.5)' C TO 'IF(XDATA(L).GT.50'IN DO 170 LOOP C JUNE 2013 IM ADDED DATA IFREQ/10*0/ C SEPTEMBER 2016 GLAHN MODIFIED COUNT FOR REL FREQ AT DO 220 C C PURPOSE C TO POSTPROCESS SKY COVER IN PERCENT AT STATIONS. C THE INTENT IS TO HAVE OPAQUE SKY (THAT CAN C BE CHECKED FOR CONSISTENCY WITH CEILING HEIGHT), BUT C SOME PREVIOUS DATA MAY BE TOTAL SKY. THIS WILL NOT C CHANGE THE COMPUTATIONS HEREIN. C C WHEN THERE IS A CEILING HEIGHT IN HUNDREDS OF FEET, C THEN IF THE SKY IS < 50 PERCENT, IT IS MODIFIED TO C A VALUE DETERMINED BY THE CLOSEST NEIGHBOR WITH C SKY COVER > 50 PERCENT. C C THE VARIABLE IS THEN SCALED BY FACTOR=CONST*10.**NSCAL. C SKYMBO IS CALLED FROM SCLSKY. C C DATA SET USE C KFILDO - UNIT NUMBER FOR OUTPUT (PRINT) FILE. (OUTPUT) C KFIL10 - UNIT NUMBER FOR INTERNAL RANDOM ACCESS STORAGE. C (INPUT) C C VARIABLES C KFILDO = UNIT NUMBER FOR OUTPUT (PRINT) FILE. (INPUT) C KFIL10 = UNIT NUMBER FOR INTERNAL RANDOM ACCESS STORAGE. C (INPUT/OUTPUT) C NDATE = DATE/TIME, YYYYMMDDHH, OF THE RUN. C (INPUT) C ID(J) = ID OF VARIABLE TO PROVIDE DATA FOR ANALYSIS C (J=1,4). (INPUT) C IDPARS(J) = THE PARSED, INDIVIDUAL COMPONENTS OF THE C VARIABLE ID'S CORRESPONDING TO ID( ,N) C (J=1,15), (N=1,ND4). C J=1--CCC (CLASS OF VARIABLE), C J=2--FFF (SUBCLASS OF VARIABLE), C J=3--B (BINARY INDICATOR), C J=4--DD (DATA SOURCE, MODEL NUMBER), C J=5--V (VERTICAL APPLICATION), C J=6--LBLBLBLB (BOTTOM OF LAYER, 0 IF ONLY C 1 LAYER), C J=7--LTLTLTLT (TOP OF LAYER), C J=8--T (TRANSFORMATION), C J=9--RR (RUN TIME OFFSET, ALWAYS + AND BACK C IN TIME), C J=10--OT (TIME APPLICATION), C J=11--OH (TIME PERIOD IN HOURS), C J=12--TAU (PROJECTION IN HOURS), C J=13--I (INTERPOLATION TYPE), C J=14--S (SMOOTHING INDICATOR), AND C J=15--G (GRID INDICATOR). C (INPUT) C JD(J) = THE BASIC INTEGER VARIABLE ID'S (J=1,4) C (N=1,ND4). C THIS IS THE SAME AS ID(J,N), EXCEPT THAT THE C PORTIONS PERTAINING TO PROCESSING ARE OMITTED: C B = IDPARS(3, ), C T = IDPARS(8,), C I = IDPARS(13, ), C S = IDPARS(14, ), C G = IDPARS(15, ), AND C THRESH( ). C NOT ACTUALLY USED. (INPUT) C IDCIG = THE FIRST WORD OF THE CEILING HEIGHT TO C USE IN CHECKING CONSISTENCY. IT COMES FROM C U405A.CN ITABLE(1,7). (INPUT) C XDATA(K) = SCALED CATEGORICAL VALUES ON INPUT; ACTUAL C VALUES IN PERCENT ON OUTPUT (K=1,NVAL). C (INPUT/OUTPUT) C XP(K) = THE X-POSITION ON THE CURRENT GRID OF STATION K C (K=1,NVAL). (INPUT. C YP(K) = THE Y-POSITION ON THE CURRENT GRID OF STATION K C (K=1,NVAL). (INPUT. C FD2(K) = WORK ARRAY (K=1,NVAL). (INTERNAL) C NVAL = NUMBER OF STATIONS OR GRIDPOINTS BEING PROCESSED. C THIS MIGHT APPLY TO A GOE GRID. THE NUMBER OF C VALUES IN XDATA( ). (INPUT) C NCAT = NUMBER OF CEILING HEIGHT CATEGORIES. (INPUT) C (NOT ACTUALLY USED.) C CONST = THE MULTIPLIER FOR SCALING. (INPUT) C NSCAL = THE POWER OF TEN FOR SCALING. (INPUT) C RAD = R(1) IN U405A, THE RADIUS OVER WHICH TO C SEARCH FOR A CLOSE STATION. (INPUT) C LSTORE(L,J) = THE ARRAY HOLDING INFORMATION ABOUT THE DATA C STORED (L=1,12) (J=1,LITEMS). (INPUT/OUTPUT) C ND9 = MAXIMUM NUMBER OF FIELDS STORED IN LSTORE( , ). C SECOND DIMENSION OF LSTORE( , ). (INPUT) C LITEMS = THE NUMBER OF ITEMS J IN LSTORE( ,L). C (INPUT/OUTPUT) C IS0(J) = MOS-2000 GRIB SECTION 0 ID'S (J=1,4). C IS1(J) = MOS-2000 GRIB SECTION 1 ID'S (J=1,21+). C IS2(J) = MOS-2000 GRIB SECTION 2 ID'S (J=1,12). C IS4(J) = MOS-2000 GRIB SECTION 4 ID'S (J=1,4). C ND7 = DIMENSION OF IS0( ), IS1( ), IS2( ), AND IS4( ). C (INPUT) C IPACK(J) = WORK ARRAY FOR GFETCH (J=1,ND5). (INTERNAL) C IWORK(J) = WORK ARRAY FOR GFETCH (J=1,ND5). (INTERNAL) C DATA(J) = WORK ARRAY FOR GFETCH (J=1,ND5). C (INTERNAL) C ND5 = DIMENSION OF IPACK( ), WORK( ), AND DATA( ). C (INPUT) C CORE(J) = SPACE ALLOCATED FOR SAVING PACKED GRIDPOINT C FIELDS (J=1,ND10). WHEN THIS SPACE IS C EXHAUSTED, SCRATCH DISK WILL BE USED. THIS IS C THE SPACE USED FOR THE MOS-2000 INTERNAL RANDOM C ACCESS SYSTEM. (INPUT) C ND10 = THE MEMORY IN WORDS ALLOCATED TO THE SAVING OF C DATA CORE( ). WHEN THIS SPACE IS EXHAUSTED, C SCRATCH DISK WILL BE USED. (INPUT) C NBLOCK = BLOCK SIZE IN WORDS OF INTERNAL MOS-2000 DISK C STORAGE. (INPUT) C NFETCH = INCREMENTED EACH TIME DATA ARE FETCHED BY C GFETCH. IT IS A RUNNING COUNT FROM THE C BEGINNING OF THE PROGRAM. THIS COUNT C IS MAINTAINED IN CASE THE USER NEEDS IT C (DIAGNOSTICS, ETC.). (OUTPUT) C L3264B = INTEGER WORD LENGTH IN BITS OF MACHINE BEING C USED (EITHER 32 OR 64). (INPUT) C ISTOP(J) = ISTOP(1) IS INCREMENTED BY 1 WHENEVER AN ERROR C OCCURS AND THE PROGRAM PROCEEDS. ISTOP(3) IS C INCREMENTED BY 1 WHEN A DATA RECORD COULD C NOT BE FOUND. (INPUT/OUTPUT) C IER = ERROR CODE. C 0 = GOOD RETURN. C OTHER VALUES FROM GFETCH. EVERY ERROR IS C FATAL FOR THIS ELEMENT. (OUTPUT) C 1 2 3 4 5 6 7 X C C NONSYSTEM SUBROUTINES USED C NONE C DIMENSION ID(4),IDPARS(15),JD(4) DIMENSION XDATA(NVAL),XP(NVAL),YP(NVAL),FD2(NVAL) DIMENSION IPACK(ND5),IWORK(ND5),DATA(ND5) DIMENSION IS0(ND7),IS1(ND7),IS2(ND7),IS4(ND7) DIMENSION LSTORE(12,ND9) DIMENSION CORE(ND10) DIMENSION FREQ(10),IFREQ(10),RFREQ(10),ISTOP(6),LD(4) DATA FREQ/01.,05.,12.,25.,37.,50.,68.,87.,93.,999./ DATA IFREQ/10*0/ C D CALL TIMPR(KFILDO,KFILDO,'START SKYMBO ') IER=0 C C CHECK CONSISTENCY WITH CEILING FORECASTS. IF THERE IS C A CEILING, THEN THE SKY MUST BE > 50 PERCENT COVERAGE. C CEILING HEIGHTS WERE WRITTEN BY CIGMBO WITH UNLIMITED C = 888. ALL VALUES LT 888 ARE LEGITIMATE CEILINGS. C THE VALUE COULD BE > 120 BECAUSE OF AUGMENTATION WITH C PERSISTENCE. C C RETRIEVE THE CATEGORICAL CEILING HEIGHT AT STATIONS. C IDCIG IS THE ID OF CEILING HEIGHT SANS MODEL NUMBER. C LD(1)=IDCIG+IDPARS(4) LD(2)=0 LD(3)=IDPARS(12) LD(4)=0 C BEFORE RETRIEVING CATEGORICAL CEILING HEIGHT, CHECK C TO SEE IF THE ID IS SPECIFIED IN THE U405A CONTROL C FILE. IF NOT, SKIP GFETCH AND GOTO 200. IF(LD(1).EQ.IDPARS(4).AND.LD(2).EQ.0.AND.LD(4).EQ.0)THEN WRITE(KFILDO,164)(LD(M1),M1=1,4) 164 FORMAT(/' NO CEILING HEIGHT ID SPECIFIED IN U405A CN.'/ 1 ' (ITABLE(J,7) IS ',3I10.9,I10,').'/ 2 ' MAY NOT BE AN ERROR. CONTINUTING WITHOUT', 3 ' MODIFICATION.') GO TO 200 ELSE CALL GFETCH(KFILDO,KFIL10,LD,7777,LSTORE,ND9,LITEMS, 1 IS0,IS1,IS2,IS4,ND7,IPACK,IWORK,DATA,NVAL, 2 NWORDS,NPACK,NDATE,NTIMES,CORE,ND10, 3 NBLOCK,NFETCH,NSLAB,MISSP,MISSS,L3264B,1,IERCIG) C IF(IERCIG.NE.0)THEN ISTOP(3)=ISTOP(3)+1 WRITE(KFILDO,165)(LD(M1),M1=1,4) 165 FORMAT(/' ****COULD NOT FIND CEILING HEIGHT FORECASTS', 1 3I10.9,I10,'. CONTINUING WITHOUT', 2 ' MODIFICATION.') GO TO 200 ENDIF C ENDIF C D WRITE(KFILDO,168)(DATA(K),K=1,NVAL) D168 FORMAT(/,' IN SKYMBO AT 176--DATA(K)',/,(15F8.2)) C RSQ=RAD*RAD C RAD IS THE FIRST RADIUS OF INFLUENCE FROM U405A.CN. C DO 180 K=1,NVAL FD2(K)=XDATA(K) C IF(XDATA(K).LT.9998.5)THEN C XDATA( ) IS SKY FORECAST IN FRACTIONS OF COVERAGE. C IF(DATA(K).LT.9998.5)THEN C DATA( ) IS CEILING FORECAST IN HUNDREDS OF FT. C CAN'T USE IT IF IT IS MISSING. C IF(DATA(K).LT.887.)THEN C UNLIMITED CEILING = 888. ALL OTHER NON-MISSING C VALUES REPRESENT A LEGITIMATE CEILING. C IF(XDATA(K).LT.50.)THEN C C MUST REPLACE THE SKY COVER. FIND THE CLOSEST C STATION WITH > 50 PERCENT COVERAGE, WITHIN C THE RANGE R(1) (THE FIRST RADIUS OF INFLUENCE), C AND SET THIS VALUE TO THAT VALUE. THIS IS TO C TRY FOR SPATIAL CONSISTENCY OF COVERAGE. C LSAVE=0 DISTMX=999999. C DO 170 L=1,NVAL C IF(L.EQ.K)GO TO 170 C DON'T USE THE STATION ITSELF. C IF(ABS(XP(K)-XP(L)).LE.RAD.AND. 1 ABS(YP(K)-YP(L)).LE.RAD)THEN C THIS CHECK IS MADE TO KEEP FROM C UNNECESSARILY CALCULATING SQUARES. C THIS TEST WILL ELIMINATE MOST STATIONS. C DISTSQ=(XP(K)-XP(L))**2+(YP(K)-YP(L))**2 C IF(DISTSQ.LE.RSQ)THEN D WRITE(KFILDO,169)K,XDATA(K),L,XDATA(L), D 1 LSAVE,DISTSQ D169 FORMAT(' AT 169--K,XDATA(K),L,XDATA(L),', D 1 'LSAVE,DISTSQ',I6,F8.2,I6,F8.2,I6,F8.1) C IF(DISTSQ.LT.DISTMX)THEN C STATION L IS CLOSER. C IF(XDATA(L).GT.50..AND. 1 XDATA(L).LT.9998.5)THEN C CLOSE STATION HAS > 50 PERCENT C COVERAGE, SO SAVE LOCATION. DISTMX=DISTSQ LSAVE=L ENDIF C ENDIF C ENDIF C ENDIF C 170 CONTINUE C IF(LSAVE.NE.0)THEN D WRITE(KFILDO,177)K,XDATA(K),LSAVE,XDATA(LSAVE) D177 FORMAT(/' AT 177--K,XDATA(K),LSAVE,XDATA(LSAVE)' D 1 ,I6,F10.3,I6,F10.3) FD2(K)=XDATA(LSAVE) ELSE D WRITE(KFILDO,178)K,XDATA(K) D178 FORMAT(/' AT 178--K,XDATA(K)',I6,F10.3) FD2(K)=51. C IF NO CLOSE STATION WITH > 50 PERCENT C COVERAGE, SET STATION TO 51 PERCENT. ENDIF C ENDIF C ENDIF C ENDIF C ENDIF C D WRITE(KFILDO,179)K,XDATA(K),FD2(K) D179 FORMAT(' AT 179--K,XDATA(K),FD2(K)',I6,2F10.2) C 180 CONTINUE C C TRANSFER THE MODIFIED VALUES TO XDATA( ). C DO 190 K=1,NVAL XDATA(K)=FD2(K) 190 CONTINUE C 200 CONTINUE C C COMPUTE THE FREQUENCIES. THE TEST IS ON LE, SO THE LISTED C VALUE IS THE HIGH END OF THE CATEGORY, INCLUSIVELY. C TOTAL=0. C DO 220 K=1,NVAL DO 218 J=1,10 C IF(XDATA(K).LE.FREQ(J))THEN IFREQ(J)=IFREQ(J)+1 TOTAL=TOTAL+1 GO TO 220 ENDIF C 218 CONTINUE C 220 CONTINUE C IF(TOTAL.EQ.0.)THEN WRITE(KFILDO,225) 225 FORMAT(/' ERROR IN SKYMBO. TOTAL = 0.') ISTOP(1)=ISTOP(1)+1 GO TO 240 ENDIF C DO 230 J=1,10 RFREQ(J)=IFREQ(J)/TOTAL 230 CONTINUE C WRITE(KFILDO,232)(FREQ(J),J=1,10) 232 FORMAT(/' LE PERCENT ',15F7.2,' TOTAL') WRITE(KFILDO,233)(IFREQ(J),J=1,10),NINT(TOTAL) 233 FORMAT(' FREQUENCIES',10I7,I9) WRITE(KFILDO,234)(RFREQ(J),J=1,10) 234 FORMAT(' REL FREQ ',10F7.2) C C SCALING NECESSARY ONLY WHEN CONST NE 1 OR NSCAL NE 0. C IF(CONST.NE.1.OR.NSCAL.NE.0)THEN C FACTOR=CONST*10.**NSCAL C WRITE(KFILDO,237)FACTOR 237 FORMAT(/' SCALING FACTOR =',F10.3,' IS BEING USED IN SKYMBO.') C DO 238 K=1,NVAL C IF(XDATA(K).LT.9998.5)THEN XDATA(K)=XDATA(K)*FACTOR ENDIF C 238 CONTINUE C ENDIF C 240 RETURN END