Miller, Marvin E.

Weaver, C. R.

KnowledgeBank at OSU

RESEARCH CIRCULAR 155 A Precipitation Probability Computer Program C.R. WEAVER MARVIN MILLER OHIO AGRICULTURAL RESEARCH AND DEVELOPMENT CENTER Wooster, Ohio NOVEMBER 1967 CONTENTS * * * I ntrod uct1on _______________________ -- - __ - - - - _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 3 Operating Procedure ___________ - ___ - ___ - - - - - - - - _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 3 Input Cards-------------------------------------------------- 3 Program Checks and Alternatives _________________________________ 4 Literature Cited _____________________________ - _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 4 Examples of Data Input_ ____________________________________________ 5 Example of Output_ _____________________ - - ____ - _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 6 Fortran Source Program ________________ - ___ - - - ___ - _ - _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 6 AGOEX 074 11-67-2.SM A Precipitation Probability Computer Program 1 C. R. WEAVER and MARVIN MILLER2 INTRODUCTION The following computer program was written and compiled in PDQ FORTRAN (1620 General Program Library 2.0.031) for the purpose of com-puting amounts of precipitation associated with vari-ous prdbability levels. It is similar to that written by McKay ( 4) but with several modifications to allow it to operate over a wider variety of values. The input data consists of a number of years of record of precipitation. This particular program operates on monthly and annual precipitation. The assumption is made that precipitation amounts fol-low the incomplete gamma distdbution (5). Two parameters of the gamma distribution, gamma and beta, are estimated from the data. For gamma values of 100 or less, precipitation corresponding to particular percentage points of the incomplete gam-ma distribution are computed. For gamma values of more than 100, precipitation corresponding to per-centage points of the normal distribution are used (5). The incomplete gamma distriibution may be ex-pressed: P(m,x) b m-1 e -b d b, where m the gamma parameter, b = X/{3 ((3 = beta parameter), and x the value of the precipitation associated with the prdbaibility P ( m, x) . Thom ( 6) showed that x for a given probaibility may be estimated by a Newton approximation as fol-lows: Letj=i-1 x· J where symbols represent values as defined above. Thom ( 6) also prepared direct and inverse tables of the gamma integral. In the computer program, the initial estimate of Xo equals the gamma parameter minus one. The approximation begins at that point. The quantity in •the brackets is defined in the program as S and is estimated to the desired accuracy by an iterative procedure. The IBM 1620, which was used for develop-ment of this program, is limited to numbers of the magnitude 10 ± 49. Gamma functions of numbers over 36 resulted in operations beyond the capability of the computer. Therefore, the terms f (m) and X m· 1 were expressed as log base e. Tables of the gamma function were used and linear interpolation between gamma, gamma + 1/a, gamma + % , and gamma + % was used if the gamma parameter fell between 36 and 100. (T1ables of log base 10 of the gamma function ( 1 ) were converted to log base e.) If the gamma parameter was less than 36 an algorithm due to Collinge (2) was used to com~ute the gamma function. OPERA TING PROCEDURE Input Cards 1. A card with 13 A3 monthly name abbrevia-tions. These abbreviations appear on the output. 2. A card with five F8.5 values of standard de-viation corresponding to the normal curve probabilities which represent the last six rainfall probability levels. 1The authors wish to express appreciation to G. F. McKay, En· vironmental Science Services Administration, National Weather Records Center, Asheville, North Carolina, for providing the Fortran IV pro-gram from which this routine was developed. •statistician, Ohio Agricultural Research and Development Center; and State Climatologist, Environmental Science Services Administration, Weather Bureou. + X 1• = X • J -m + x + m+l (m+l) (m+2) (m+l) (m+2) (m+3) Pr(m)e X· m-1 J 3 3. Sixty-five cards containing X-35, the value of log base e of the gamma functions of X = 36 to 100, X + 0.33, X + 0.5, and X + 0.67 (format is I3,4El 6.8). 4. For each year of record: a) 12 cards with station, year, month, and monthly precipi-tation; and b) one card with station, year, and annual precipitation. The format is cc 1-6 = Station; 8-9 =Year; 10-11 =Month; 43-46 = Total precipi-tation on monthly cards; 65-68 = Total precipitation on anuual cards; 80 = 3 on annual cards. 5. Following the data one card with -1 in cc 79-80. 6. During the input, the following checks are performed on the tables and the data. a. X-35 (X = 36 to 100) is checked for sequence. A sequence error re-quires that the entire table be re-en-tered after sequence correction. b. During the data reading, the station num!ber of each monthly and annual card is checked against the first month. No year to year check of station is made. c. A check of the sequence of the monthly cards is also made. Any error in b or c sends the program back to read the entire year which is in error. The input data cards being used in this program are produced by two weather summary programs developed by Lytle and Martens (3). The FOR-TRAN program could easily be modified to utilize other input formats. Examples of the inputs are listed on page 5. Program Checks and Alternatives 1. Two iterative loops may cause difficulty. The first is connected with the hypergeometric series which is used to evalurute S. The change in the esti-mate of S is called T in the program and io; tested 4 against TEST. The value of TEST is arbitrarily set at 0.0001. The limit for the number of iterations in computing S is 1000. In the data processed to date, 10 to 30 iterations have been sufficient. SWITCH 1 on will print out current values for the month, the probability level index, the number of iterations, the estimate of S, the estimate of T, and TEST. SWITCH 2 on allows a manual entry on the typewriter of a new value of TEST. An Fl0.0 for-mat is used so one must be careful to type in 10 characters to assure replacement of other digits which may be in memory. 2. The second iterative estimate is performed on rainfall divided by the beta parameter. This esti-mate is called X1 in the program. It is possible for this iteration to proceed indefinitely. After 10 itera-tions, a message to that effect is automatically typed out. SWITCH 3 on will print out current values for the month, the proba!bility level index, the number of iterations, the estimate of X1, the value of the change in X1, and DEST which is the criterion against which DX (change in X 1 ) is tested. SWITCH 4 on allows entry of a new value of DX. In approximately 30 trials of actual Ohio wea-ther data, the program never exceeded 10 iterations. LITERATURE CITED 1. Abramowitz, Milton and Irene A. Stegun. 1965. Handbook 1of mathematical functions with formu-las, graphs, and mathematical tables. National Bureau of Standards, Applied Mathematics Series 55. 2. Collinge, Robert M. 1961. Gamma function. Asso-ciatiion of Computing Machinery, Algorithm 31. 3. Lytle, W. F. and Harry A. Martens. 1964. SPS weather summary SDSU. Personal communica-tion. 5. Thom, H. C. S. 1958. A note ion the gamma dis· tribution, Monthly Weather Review 86 (4): 117-122. 6. Thom, H. C. S. 1967 (in press). Direct and inverse tables of the gamma distribution. ESSA-EDS Tech. Memo No. 6. EXAMPLES OF DATA INPUT Monthly Name Card JANFEB~ARAPRMAYJUNJULAUGSEPUCTNOVDECANN Normal Curve Parameter Card 0.25335 o.52440 o.84162 1.28155 1.64485 Log Base e of Gamma Functions X-35 GAMMA(Xl GAMMA(X+l/3) GAMMA(X+l/2) 1 .92136373E 02 .93327790E 02 .93924663E 02 2 .95719898E 02 e96920533E 02 .97521983[ 02 3 .99330824E 02 el0054043E 03 .10ll4633E 03 4 .1Q296842E 03 .1Q418675E 03 .I0479700E 03 • • • • • • • • • • • • 65 .35913496E 03 .36066892E 03 .36143631E 03 Data Follows Example of 1 Year's Data GAMMA(X+2/3) .94522307E 02 .98124184E 02 .10175297£ 03 .10540795£ 03 • • • • .36220397E 03 This Is a Test Deck and Does Not Contain All of the lnformatiion That Is Outpwt by Lytle's Program. 3312?4 3101 331234 310? 3·:q234 310? 331234 3104 331234 3105 331234 3106 331234 3107 331234 3108 3'H234 3109 ?31234 3110 3'31/?4 ?111 31>1234 3112 331234 31 Several Years of Data Would Follow. 0102 0232 0213 0543 0320 0376 0427 0488 0452 0211 0277 0375 A Trailer Card with a -1 in cc 79-80 is Required After the Data Set. 5 4016 3 EXAMPLE OF OUTPUT PRECIPITATION WITH PROBABILITY EQUAL OR LE~S THAN STATION • 0 I) .10 .20 .30 .40 .so .60 • 70 .so e90 333029 JAN 1.14 1.48 2. 00 2.44 2.86 3.30 3.78 4.35 5.08 6.22 333029 FEB .86 1.18 1.66 2.09 2.52 2.97 3.46 4.05 4.82 6.03 333029 MAR 1.56 1.95 2. 50 2.96 3.39 3.83 4.32 4.87 5.58 6.67 333029 APR 1.40 1.75 2.21 2.69 3.10 3. 52 3.97 4.49 5.16 6.19 333029 MAY 1.68 2.00 2.45 2.81 3.15 3.49 3.85 4.26 4.78 5.57 3330?9 JUN 1.52 1.87 2.35 2.75 3.13 3.52 3.93 4. 41 5.02 5.94 333029 JUL 1.76 2.16 2.71 3.17 3.60 4.03 4.50 5.04 5.73 6e77 333029 .A.UG 1.34 1.67 2.16 2.56 2.95 3.34 3.77 4e26 4.89 5.86 3'33029 SEP 1.13 1.41 1.81 2.15 2.47 2.79 3.14 3.55 4.07 4086 333029 OCT .39 .58 .89 1.18 1.47 1.79 2.15 2. 58 3.16 4e08 333029 NOV .83 1.10 1.51 1.87 2.21 2.57 2.97 3. 4'+ 4e05 5.01 333029 DEC .96 1.22 1. 61 1.93 2.24 2.56 2.91 3.31 3.83 4.63 333029 ANN 33.14 34.75 36.77 38.28 39.59 40.85 42.14 43. 54 45.22 47.63 FORTRAN SOURCE PROGRAM C PRECIPITATION PROBABILITIES DIMENSION GF(65) tGF3(65),GF5(65),GF6(6~),SSSC13),CN(lll CIMENSION PL(llJ,U( 35),R(llhPC(l3l DIMENSION SS(l3),SL(l3),AVf13l,AL(l3),ALA(l3),AM<l3l C SETS PROBABILITY LEVELS PL(ll=0.05 PL(lll=0.95 DO 72 J=2,10 PL(JJ=J-1 72 PL!Jl=PL(Jl/10. C SETS GAMMA FUNCTION COEFFICIENTS A=.16063118E-02 V=.51589951E-02 C=.44511400E-02 D=.72110256E-Ol E=.82111740E-Ol F=.41177419 G=.42278746 H=.99999997 C GAMMA FU~CTION TABLE INTERPOLATION POINTS Zl=.16666667 Z3=.33333333 Z5=.5 Z6=.66666667 C MONTH IDENTIFICATION .95 1.21 7.16 7.66 7.14 6.28 6.79 1.12 6.75 5.59 4~96 5e89 5.37 49.69 READ 61Q,AM(l),AM<2J,AM(3ltAMC41,AM(5),AM(6),AM(7l ,AMC8J,AM(9), lAM<lO> ,AM(ll) ,Alv\<12) ,AM(l3) 610 FORMAT(l3A3) 6 C NORMAL CURVE PARAMETERS READ 615,CN(7),(N(8) ,CN(9),(N(l0l,CN!lll 615 FORMATC5F8.5l DO 205 J=l,5 K=J2-J 205 CN(JJ=CNCKl*C-1.0l CN!6l=o.o C READS GAMMA FUNCTION TABLES 208 DO 206 J=l,65 READ 616,L,GFCJJ,Gf3(J) ,GF5(JJ,GF6(J) 616 FORMAT!I3,4El6.8l IFCL-Jl97t206,97 97 PRINT 617 617 FORMATC17HGF SEQUENCE ERROR) PAUSE GO TO 208 206 CONTINUE 1 DO 2 J=1'13 SSCJ>=o.o SSS(Jl=o.o 2 SLCJl=O·O CT=o.o C READ 12 MONTHLY AND ONE ANNUAL DATA C INPUT ERROR MESSAGE REQUIRES READING ALL CARDS FOR YEAR IN C ERROR AGAIN AFTER CORRECTING SEQUENCE 9 DO 7 J=l,12 READ 601,STA,STB,IY,IM,PCCJ),IB 601 FORMATC2A3,1x,2I2,31x,F4.2,32X,I2) IFCJ-1)3,3,4 3 IFCIBJ5'6'6 C TRAILER WITH -1 IN CC 79-80 CAUSES BREAK TO STATEMcNT ~ 6 STAS=STA STBS=STB IYS=IY GO TO 7 4 IFCIM-JJ99,8,99 99 PRINT 602 602 FORMATC9HMONTH ERRl PAUSE GO TO 9 8 IFCSTAS-STAl98,6Q,98 60 IFCSTBS-STB)98,61,98 61 IF!IYS-IYl98,7,9a 98 PRINT 603 603 FORMATCllHSTATION ERR) PAUSE GO TO 9 7 CONTINUE READ 604,STA,STB,JY,PC(l3l 604 FORMATC2A3,1x,12,55x,F4.2) IFCSTAS-STAl98,62,98 62 IFCSTBS-STBl98,63,98 63 IFCIYS-IYl98,64,98 7 C ADDS SQUARES AND LOGS 64 DO 65 J=ltl3 SS(Jl=SS(J)+PC(J) SSS(Jl=SSS!J)+(PC(Jl*PC(J)l IFCPC(Jll82t82t65 82 PC<J>=OeOOl 65 SL(JJ=SL<Jl+LOGF(PC(J)J CT=CT+leO GO TO 9 C PUNCHES HEADERS 5 PUNCH 608 608 FORMAT(l8Xt49HPRECIPITATION WITH PROBABILITY EQUAL OR LESS THAN, 1 13H l PUNCH 609tPLClltPLC2ltPLC3J,PLC4ltPL(5ltPLC6ltPL(7ltPLC8),PL(9), 1 PL(lOJ,pL(ll) 609 FORMAT(7HSTATIONt3XtllF6.2l PUNCH 605 C THIS LOOP DOES 12 MONTHLY AND AN ANNUAL SET OF PROBABILITIE~ BEGIN TRACE DO 66 J=ltl3 AV(JJ=SS<Jl/CT AL(Jl=LOGFCAVCJll ALACJl=SL(J)/CT Z=l2•0*CALCJl-ALACJJ l C CALCULATES GAMMA PARAMETER GG=(6.Q+SQRTFC36.0+C4.Q*Zl ll/(2•0*Zl C USES NORMAL CURVE IF GAMMA OVER 100 lF(GG-lOQ.)200t200t201 201 SSSCJl=SSSCJl-CAV(Jl*SS(Jll SD=SQRTFCSSSCJl/(CT-laOll DO 203 I=l,11 RCI>=AV{J)+CCNCil*SD) IFCRCill202,203t203 202 R(ll=O·O 203 CONTINUE GO TO 207 200 X=AV(J) B=X/GG P=GG-laO L=P+l.O IFCL-ll13,13tl0 10 UC 1> =GG C INTERPOLATES GAMMA FUNCTION OF GAMMA FOR GG BETWEEN 36 AND 100 IF!GG-36ell09,110tll0 110 K=GG-35. SB=K DEC=GG-SB-35. IF!DEC-Z6llll,113tll3 113 M=K+l GAM=GF6CKl+CC <DEC -Z6l/Z3J*CGF!MJ-Gf6(Kll l GO TO 103 111 IF<DEC-Z5)114tll6tll6 116 GAM=GF5!Kl+Cf CDEC•Z5l/Zll*CGF6CKJ-GF5(Kll > GO TO 103 114 IFCDEC-Z3lll7,120tl20 8 120 GAM=GF3CKl+C(CDEC-Z3J/Zll*{GF5(Kl-GF3(K)l l GO TO 103 117 GAM=GFCKJ+(CDEC/Z3l*(GF3(Kl-GF(Kll l GO TO 103 C NUMERICAL APPROXIMATION OF GAMMA FUNCTION OF GAMMA FOR GG LESS THAN 36 109 DO 11 K=ltL KK=K+l 11 UCKKl=UCKl-1.0 PD=l·O DO 12 K=2tL 12 PD=PD*UCKl Y=U ( Ll PD=PD/Y GO TO 14 13 Y=GG PD=l.O/Y 14 GAM=PD*( ( (( ((((A*Y+Vl*Y+Cl*Y+Dl*Y+El*Y+Fl*Y+Gl*Y+Hl/CY+l.Ol I C ESTIMATES THE VALUE OF Xl WHICH MULTIPLIED BY BETA GIVES C RAINFALL PROBABLITY 10 3 DO 2 0 I= 1'11 NEST=O TEST=O•OOOl DEST=0.0004 MEST=lO XO=P IFCGG-l.Ol5ltlOOtlOO 100 IFCXOl50,50tl5 51 X0=0.0001 GO TO 15 50 XO=O.Ol 15 Z=2•0 NEST=NEST+l T=XO/CP+Zl S=T+l•O C ITERATIVE ESTIMATE OF S c 74 611 c 76 612 75 16 17 101 122 DO 16 L=l,1000 Z=Z+leO T=<T*XO)/(P+Zl ~WITCH 1 INTERROGATES VALUE OF T IFCSENSE SWITCH 1)74,75 PRINT 6lltJtitLtStTtTEST FORMATC3I5t3El6.8l PAUSE SWITCH 2 ALLOWS ENTRY OF NEWT CRITERION CTEST> IFCSENSE SWITCH 2l76t75 ACCEPT 612tTEST FORMAT(FlO.Ol S=S+T IFCT-TESTll7tl7tl6 CONTINUE DL=PL(l) IFCDL>41tl0ltl01 IFCGG-36.)12ltl22tl22 ZUD=EXPFCGAM-CP*LOGf(XOlll GO TO 105 9 121 ZUD=EXPF(LOGF(GAMJ-(P*LOGF(XO) l) 105 Xl=XO- fXO/GGl*S +(DL*ZUD*EXPF(XOl l DX=Xl-XO C SWITCH 3 INTERROGATES VALUE OF XO IF<SENSE SWITCH 3)77,78 77 PRINT 611,J,I,NEST,Xl,DX,DEST PAUSE C SWITCH 4 ALLOWS ENTRY OF NEW CRITERION FOR DX (DEST) IFCSENSE SWITCH 4l79,78 79 ACCEPT 612,DEST 78 XO=Xl IF<XOl41,18,19 19 IF(ARSF<DX>-DESTll8,18,80 80 IF<NEST-MEST)l5,81,81 C PRINTS WARNING EVERY TEN XI ITERATIONS SO THAT OPERATOR CAN INTERCEDE 81 PRINT 614,J,I,NEST 614 FORMAT<2I5t17H XI ITERATIONS=,I5l MEST=MEST+lO GO TO 15 41 Xl=O.Q 18 R(Il=Xl*B 20 CONTINUE END TRACE 207 IFiJ-13l7Q,71,71 71 PUNCH 605 605 FOR MA TC BOX) 70 PUNCH 606,STAS,STBS,AM(J),R(~J,R(2ltR(3),R(4),R(5l,R(6ltR(7J,R(8), 1 R(9),R(l0J,R(ll) 606 FORMAT<2A3,1XtA3,llF6.2J 66 CONTINUE PUNCH 607 607 FORMATC79Xt1H-l GO TO 1 END 10 7~ State 'l6 ~ ea~ lo't ;'/~ ~e4ealid ad Z)~ Ohio's major soil types and climatic conditions are represented at the Re-search Center's 11 locations. Thus, Cen-ter scientists can make Held tests under conditions similar to those encountered by Ohio farmers. Research is conducted by 13 depart-ments on more than 6000 acres at Center headquarters in Wooster, nine branches, and The Ohio State University. Center Headquarters, Wooster, Wayne County: 1918 acres Eastern Ohio Resource Development Cen-ter, Caldwell, Noble County: 2053 acres Mahoning County Experiment Farm, Can-field: 275 acres Muck Crops Branch, Willard, Huron County: 15 acres North Central Branch, Vickery, Erie Coun-ty: 335 acres Northwestern Branch, Hoytville, Wood County: 247- acres Southeastern Branch, Carpenter, Meigs County: 330 acres Southern Branch, Ripley, Brown County: 275 acres Vegetable Crops Branch, Mari·etta, Wash-ington County: 20 acres Western Branch, South Charleston, Clark County: 428 acres 

A precipitation probability computer program

https://kb.osu.edu/bitstream/1811/70605/1/OARDC_research_circular_n155.pdf

A precipitation probability computer program

Abstract

Similar works

Full text

Available Versions

KnowledgeBank at OSU