C wlf_GEM6.f C ################################################################## C ################################################################## C ###### ###### C ###### ###### C ###### PROGRAM FOR DAILY WEATHER SIMULATION ###### C ###### IN THE CONTIGOUS UNITED STATES ###### C ###### ###### C ###### ARS - 114 ###### C ###### (c) July 1994 ###### c ##### Feburay 2000 ###### C ###### ###### C ###### UNITED STATES DEPARTMENT OF AGRICULTUE ###### C ###### AGRICULTURAL RESEARCH SERVICE. ###### C ###### All rights reserved. ###### C ###### ###### C ################################################################## C ################################################################## PROGRAM GEM6 C####################################################################### C C PURPOSE: C C This program provides precipitation probabilities and simulates C data for daily precipitation, maximum and minimum temperature, and C solar radiation for an n-year period at a given location within the C contiguous United States. Wind Speed and Dewpoint were added in c Feb. 2000. The model is designed to preserve the dependence in c time, the internal correlation, and the seasonal characteristic c that exist in actual weather data. Daily maximum and minimum c temperture, dew point, wind speed, and solar radiation data are c simulated using a weakly stationary generating process conditioned c on the precipitation process described by a Markov chain-mixed c exponential model. c C Paratmeters for special stations within a region can be accessed C directly. An estimation routine(s) will be added for estimation of c data (and / or parameters) for points between existing stations. c C The seasonal variations of parameters are described by the Fourier C series. Information on the type of equipment needed to run the model C and an example of running the model are provided. C C####################################################################### C C AUTHOR: C.L.Hanson,K.A.Cumming,D.A.Woolhiser and C.W.Richardson C C Modify history: C Yunyun Lu (November,1994) C transfer program from QBASIC to FORTRAN 77 C C Y. Lu (Dec. 1994) C added the minimal documentation C C###################################################################### C Will Frymire (Dec 1999) c transefered to PC version: c No graphics. c Input Name of Station location file (i.e. "Boise" for Boise Idaho) c Added comments c Output file 'filename.out' starts in January.(i.e. Boise.out) C####################################################################### C C Variable Declarations. C C####################################################################### C C OUTPUT: C C p00 Transition probability of a dry day following a dry day C p10 Transition probability of a dry day following a wet day C beta The mean of the smaller exponential distributions C delta The mean of the larger exponential distributions C albar The weighting function C C anp Expected annual precipitation (inches) C sum1 number of wet days C C p Simulated M years of daily rainfall C tmax Simulated M years of daily max. temperature C tmin Simulated M years of daily min. temperature C rad Simulated M years of daily total radiation C C####################################################################### C C Variable Declarations. C C####################################################################### C implicit none integer idyr,m parameter (idyr=365,m=360) common /const1/pi common /sited2/delta(idyr),p00(idyr),p10(idyr),beta(idyr),albar common /rainf1/sum1,anp common /tem11/tmad(idyr),txs(idyr),txm1(idyr),txs1(idyr), !dmad(idyr),dxs(idyr),dxm1(idyr),dxs1(idyr), !wmad(idyr),wxs(idyr),wxm1(idyr),wxs1(idyr), !tnm(idyr),tns(idyr),tnm1(idyr),tns1(idyr), !rm0(idyr),rs0(idyr),rm1(idyr),rs1(idyr), !rm(idyr),chipre(6) common /temp2/p,tmax,tmin,dew,wind,rad common /temp3/fran(12) common /temp4/a(6,6,12),b(6,6,12) real albar real delta,p00,p10,beta real pi real p,tmax,tmin,dew,wind,rad real sum1,anp real tmad,txs,txm1,txs1,tnm,tns,tnm1,tns1 real rm0,rs0,rm1,rs1,rm,chipre real dmad,dxs,dxm1,dxs1 real wmad,wxs,wxm1,wxs1 real fran real a,b C C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C C Beginning of executable code... C C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C write(6,'(/ 22(/5x,a)//)') :'###############################################################', :'###############################################################', :'##### #####', :'##### Welcome to #####', :'##### #####', :'##### UNITED STATES WEATHER SIMULATOR GEM6 #####', :'##### #####', :'##### ARS - 114 #####', :'##### (c) July 1994 #####', :'##### FEB. 2000 #####', :'##### #####', :'##### This program provides easy access to informate #####', :'##### daily precipitation, maximum and minimum #####', :'##### temperture, and solar radiation. #####', :'##### #####', :'##### UNITED STATES DEPARTMENT OF AGRICULTUE #####', :'##### AGRICULTURAL RESEARCH SERVICE. #####', :'##### All rights reserved. #####', :'##### #####', :'###############################################################', :'###############################################################' pi=atan(1.)*4 CALL INPUT_PARAMETERS ! input fourier coefficients for c ! each station END C ################################################################## C ################################################################## C ###### ###### C ###### SUBROUTINE INPUT_PARAMETERS ###### C ###### ###### C ###### (c) 1994 ###### C ###### ###### C ################################################################## C ################################################################## SUBROUTINE INPUT_PARAMETERS C C####################################################################### C C PURPOSE: C C Input Fourier coefficients for each station. C C Parameter interpolation for user chosen station. One option used. C (Choose Closest station) C C Other options to be implemented at a later date: (after Feb. 2000) C (1) to use arithmetic average of parameters from stations C within a radius of 100 miles. C (2) to use the parameters of the nearest neighbor. C C####################################################################### C C AUTHOR: C.L.Hanson,K.A.Cumming,D.A.Woolhiser and C.W.Richardson C C Modify history: C Yunyun Lu (November,1994) C transfer program from QBASIC to FORTRAN 77 C set up the algorism for distance oon curved earth surface C C Y. Lu (Dec. 1994) C added the full documentation C C####################################################################### C C INPUT: (Inputs from AGUA46 ) This version of GEM uses only 3 harmonics C C Total 360 stations in US region C C x station's longitude. C y station's latitude. C stname station's name. C z parameter (p00,p10,beta,mu) matrix for C MCME model (upto 6 harmonics). C 1. Annual mean C 2. Amplitude of first harmonic C 3. Phase angle of first harmonic C 4. Amplitude of second harmonic C 5. Phase angle of second harmonic C 6. Amplitude of third harmonic C 7. Phase angle of third harmonic C 8. Amplitude of fourth harmonic C 9. Phase angle of fourth harmonic C 10. Amplitude of fifth harmonic C 11. Phase angle of fifth harmonic C 12. Amplitude of sixth harmonic C 13. Phase angle of sixth harmonic C zz weighting function C C OUTPUT: C C average parameter for the station C C p00 Transition probability of a dry day following a dry day C p10 Transition probability of a dry day following a wet day C beta The mean of the smaller exponential distributions C delta The mean of the larger exponential distributions C albar The weighting function C C WORKING ARRAY: C C th C za C r distance between center and stations C C####################################################################### C C Variable Declarations. C C n nh (Aqual coefficients array size (n by,nh) ) C C####################################################################### C parameter(m=360,n=4,nh=7,idyr=365) !(nh-1)/2= # harmonics real th(idyr),za(n,nh) ! nh= 13 for 6 harmonics character stname(m)*20, stna(3)*20,char1 common /const1/pi common /sited2/delta(idyr),p00(idyr),p10(idyr),beta(idyr),albar common /rainf1/sum1,anp common /temp3/fran(12) common /temp4/a(6,6,12),b(6,6,12) character*25 site,temp,matrix,gem6out,filename, AguaHeader common/filname/site,temp,matrix,gem6out,filename c####################################################################### C C VarName filename.extension PROGRAM Comments C site Location.site - AGUA generated precipitation file C temp " .temp - GENPAR generated annual means (coefficients) C matrix " .max - GENPAR generated A & B matrix (March 1) C gem6out " .out - GEM6 X num years of generated weather (Jan1) C " .txt - SAMSON (CD) --> SAMFORMT--> *.txt file C filename - location name of file (i.e. Boise for Boise.txt) C *.site contains coeficients for stations c####################################################################### c Read data from file "sites.dat" by station name,latitude, c longitude and parameter matrix. Calculate the distance between c center (provided by user) and each station, store it in array c r(m). c####################################################################### DO 10 i=1,idyr th(i)=0. 10 CONTINUE albar=0. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc write (6,*)'Enter Location file name (no extension-no blanks) ' read(5,*) filename c c Fortran pads the character string 'filename' with blanks. Check for c the 1st blank character. This is the end of the file name string. Look c for 'filename.site' etc. c ilen = index(filename,' ') ! string position of 1st ' ' c site=filename(1:ilen-1)// '.site' ! from AGUA46 - ppt output c temp = filename(1:ilen-1) // '.temp' ! from GENPAR - annual means c matrix = filename(1:ilen-1) // '.max' ! GENPAR - A&B matrix c gem6out = filename(1:ilen-1) // '.out' ! GEN6 output file c c####################################################################### c Begining the reading process c####################################################################### ccc open(2,file='boise61a.site',status='old') ! Input SITES (Agua) c open(2,file=site,status='old') ! Input SITES (Agua) ccc plat1=43.+43./60. !HARD CODED for BOISE c read(2,*)AguaHeader read(2,*)plat1 read(2,*) ((za(i1,i2),i1=1,n),i2=1,nh) ! Agua coefficients read(2,*) albar ! read weighting coefficient close(2) c c####################################################################### c Calculate the seasonal variations (see orange book Eq. (3) Pg 2 c####################################################################### c t=2*pi/float(idyr) DO 23 i1=1,n DO 24 i=1,idyr,4 th(i)=za(i1,1) IF(za(i1,2).gt.0.0001) ! th(i)=th(i)+za(i1,2)*sin(t*i+za(i1,3)) IF(za(i1,4).gt.0.0001) ! th(i)=th(i)+za(i1,4)*sin(t*i*2+za(i1,5)) IF(za(i1,6).gt.0.0001) ! th(i)=th(i)+za(i1,6)*sin(t*i*3+za(i1,7)) cccCCCCCCCCCCCCC for 6 harmonics CCCCCCCCCCCCCCcc ccc IF(za(i1,8).gt.0.0001) ccc ! th(i)=th(i)+za(i1,8)*sin(t*i+za(i1,9)) ccc IF(za(i1,10).gt.0.0001) ccc ! th(i)=th(i)+za(i1,10)*sin(t*i*2+za(i1,11)) ccc IF(za(i1,12).gt.0.0001) ccc ! th(i)=th(i)+za(i1,12)*sin(t*i*3+za(i1,13)) c ccccccccccccccccccccccccccccccccccccccccccccc IF(i.gt.1) THEN DO 27 k=1,3 th(i-k)=th(i-4)+(th(i)-th(i-4))*(4-k)/4. 27 CONTINUE END IF 24 CONTINUE DO 26 i=1,idyr IF(i1.eq.1) p00(i)=th(i) IF(i1.eq.2) p10(i)=th(i) IF(i1.eq.3) beta(i)=th(i) IF(i1.eq.4) delta(i)=(th(i)-albar*beta(i))/(1.-albar) ! the calculation of delta is according to orange ! book Eq. (6), page 3 26 CONTINUE 23 CONTINUE write(6,'(/ 22(/5x,a)//)') :' ', :' ', :'Do you want to calculate average annual precipitation (Y/N)? ' read(5,'(a1)') char1 IF(char1.eq.'N'.or.char1.eq.'n') THEN write(6,'(5x,a)') 'hit "return"/"enter" key please.' read(5,*) GOTO 30 END IF CALL ANNUAL_RAINFALL ! average annual rainfall write(6,'(/ 22(/5x,a)//)') :'If you know average annual precipitation please enter: ', :'(inches) in form of xxx.xx. ' read(5,'(f10.5)') avp c read(5,*) avp IF(avp.eq.0.0) GOTO 30 c####################################################################### c A NEWTON-RAPHSON iterative procedure is used to adjust c alfarbar and P10 bar, so that the theoretical mean (from c interpolation isohyetal map) within +/- 0.1 percent of c the known average annual precipitation c####################################################################### 2 eer=avp-anp IF(abs(eer).lt.0.001) GOTO 30 dedp1=(1-za(1,1))/(1+za(2,1)-za(1,1))**2 ! *(albar*za(3,1)+(1-albar)*za(4,1))*idyr ! orange book Eq. (25), page 11 deda=-(1-za(1,1))*(za(3,1)-za(4,1))*idyr/(1+za(2,1) ! -za(1,1)) ! orange book Eq. (26), page 12 da=-0.33*eer/deda ! orange book Eq. (23), page 11 dp1=-0.33*eer/dedp1 ! orange book Eq. (24), page 11 za(2,1)=za(2,1)+dp1 albar=albar+da DO i=1,idyr p10(i)=p10(i)+dp1 END DO CALL ANNUAL_RAINFALL GOTO 2 30 CONTINUE write(6,'(/ 22(/5x,a)//)') :' ', :' ', :'Do you want to estimate the probability of x or less inches of', :'precipitation over N days ? (Y/N) ' read(5,'(a1)') char1 IF(char1.eq.'Y'.or.char1.eq.'y') THEN CALL TOTAL_RAINFALL ! for n-day total precipitation ELSE write(6,'(5x,a)') 'hit "return"/"enter" key please.' read(5,*) END IF write(6,'(5x,a)') :'Do you want to start the simulation process? (Y/N) ' read(5,'(a1)') char1 IF(char1.eq.'Y'.or.char1.eq.'y') THEN C CC write(6,'(5x,a)') CC :'Simulated data will be saved in file "Yourfilename.out" ' C CALL SIMULATION(plat1) ! for simulation of p, tmax, tmin and radiation END IF 31 RETURN END c ################################################################## c ################################################################## c ###### ###### c ###### SUBROUTINE ANNUAL_RAINFALL ###### c ###### ###### c ###### (c) 1994 ###### c ###### ###### c ################################################################## c ################################################################## SUBROUTINE ANNUAL_RAINFALL c c####################################################################### c c PURPOSE: c c Calculate average annual (365 days) rainfall c c####################################################################### c c AUTHOR: C.L.Hanson,K.A.Cumming,D.A.Woolhiser and C.W.Richardson c c Modify history: c Yunyun Lu (November,1994) c transfer program from QBASIC to FORTRAN 77 c c Y. Lu (Dec. 1994) c added the full documentation c c####################################################################### c c INPUT: c c p00 Transition probability of a dry day following a dry day c p10 Transition probability of a dry day following a wet day c beta The mean of the smaller exponential distributions c delta The mean of the larger exponential distributions c albar The weighting function c c OUTPUT: c c anp Expected annual precipitation (inches) c sum1 Number of wet days c c WORKING ARRAY: c c pwet unconditional probability of day n being wet. c c####################################################################### c c Variable Declarations. c c####################################################################### c parameter(n=4,nh=13,idyr=365) common /const1/pi common /sited2/delta(idyr),p00(idyr),p10(idyr),beta(idyr),albar common /rainf1/sum1,anp character*25 site,temp,matrix,gem6out,filename common/filname/ site,temp,matrix,gem6out,filename C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C C Beginning of executable code... C C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% pwet=(1-p00(idyr))/(1+p10(idyr)-p00(idyr)) ! orange book Eq. (4), page 2 sum1=0. sum2=0. DO 30 i=1,idyr pwet=pwet*(1-p10(i))+(1-pwet)*(1-p00(i)) sum1=sum1+pwet sum2=sum2+pwet*(albar*beta(i)+(1-albar)*delta(i)) 30 CONTINUE anp=sum2+sum1*0.01 print*,'Expected annual precipitation', anp,'inches' print*,'Number of wet days=',sum1 RETURN END c ################################################################## c ################################################################## c ###### ###### c ###### SUBROUTINE SIMULATION ###### c ###### ###### c ###### (c) 1994 ###### c ###### ###### c ################################################################## c ################################################################## SUBROUTINE SIMULATION(PLAT1) c####################################################################### c c PURPOSE: c c Simulate M years of daily rainfall data beginning Mar.1 c c####################################################################### c c AUTHOR: C.L.Hanson,K.A.Cumming,D.A.Woolhiser and C.W.Richardson c c Modify history: c Yunyun Lu (November,1994) c transfer program from QBASIC to FORTRAN 77 c c Y. Lu (Dec. 1994) c added the full documentation c c Greg Johnson (Dec. 12, 1994) c added the option to correct the simulated temperature in these c two cases: 1) tmax(day i) < tmin(day i-1); c 2) tmin(day i) > tmax(day i-1); c c to 1) tmax(day i) = tmin(day i-1); c 2) tmin(day i) = tmax(day i-1). c####################################################################### c c INPUT: c c p00 Transition probability of a dry day following a dry day c p10 Transition probability of a dry day following a wet day c beta The mean of the smaller exponential distributions c delta The mean of the larger exponential distributions c albar The weighting function c c OUTPUT: c c p Simulated M years of daily rainfall c tmax Simulated M years of daily max. temperature c tmin Simulated M years of daily min. temperature c rad Simulated M years of daily total radiation c c WORKING ARRAY c c tmad mean of tmax -- dry days c txs standard deviation of tmax -- dry days c txm1 mean of tmax -- wet days c txs1 standard deviation of tmax -- wet days c tnm mean of tmin -- wet and dry days c tns standard deviation of tmin -- wet and dry days c rm0 mean of radiation -- dry days c rs0 standard deviation of radiation -- dry days c rm1 mean of radiation -- wet days c rs1 standard deviation of radiation -- wet days c rm max. radiation along a given latitude c chipre last step of chi(n), = chi(n-1) c chi a vector having elements that are the standardized residuals c tmax, tmin, and radiation. c c####################################################################### c c Variable Declarations. c c####################################################################### c parameter(idyr=365,ii=360) ! ii=number of stations common /const1/pi common /sited2/delta(idyr),p00(idyr),p10(idyr),beta(idyr),albar common /tem11/tmad(idyr),txs(idyr),txm1(idyr),txs1(idyr), !dmad(idyr),dxs(idyr),dxm1(idyr),dxs1(idyr), !wmad(idyr),wxs(idyr),wxm1(idyr),wxs1(idyr), !tnm(idyr),tns(idyr),tnm1(idyr),tns1(idyr), !rm0(idyr),rs0(idyr),rm1(idyr),rs1(idyr), !rm(idyr),chipre(6) common /temp2/p,tmax,tmin,dew,wind,rad common /temp3/fran(12) common /temp4/a(6,6,12),b(6,6,12) integer notdone,ycount,sw,mday(12) character stname(ii)*20,char1*4,char2*4 data mday/1,32,62,93,123,154,185,215,246,276,307,338/ character*4 mname(12) data mname/'MAR','APR','MAY','JUN','JUL','AUG', !'SEP','OCT','NOV','DEC','JAN','FEB'/ ! start from March 1 !!! c character*25 site,temp,matrix,gem6out,filename common/filname/site,temp,matrix,gem6out,filename C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C C Beginning of executable code... C C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% write(6,'(/ (/5x,a)//)') :'How many years do you want to simulate' read(5,'(i5)') m m=m+1 !simulate additional year write(6,'(/ (/5x,a)//)') :'Please give the seed for random process (-32768 -- +32768).' read(5,'(i8)') iseed leap=0 ycount=0 write(6,'(/ (/5x,a)//)') :'Do you want leap years (Y/N) ?' read(5,'(a1)') char1 IF(char1.eq.'Y'.or.char1.eq.'y') THEN write(6,'(/ (/5x,a)//)') :'Which of the first 4 years of simulation will be a leap year?' read(5,'(I2)') leap ELSE write(6,'(5x,a)') 'hit "return"/"enter" key please.' read(5,*) END IF ut=(1-p00(idyr))/(1+p10(idyr)-p00(idyr)) ! unconditional probability of day n being ! wet. (orange book,Eq.(4), Page 4) u=ran3(iseed) ! random number nw=0 ! preceding day is dry IF(u.lt.ut) nw=1 ! proceding day is wet, if u is less than ut c open(unit=21,file='x1.out',status='replace') CALL READ_123IMATION(plat1) c ccc open(unit=22,file='gem6wx.out',status='replace') !GEM6 output C (example: gem6out ='Boise.out') open(unit=22,file=gem6out,status='replace') !GEM6 output CCCCC CCCCC need to add header information CCCCC write(22,'(/ (1x,a14,6x,a15,7x,a16,6x,a3))') :'YR MO DAY prec','tmax tmin','dewpoint wind','rad' c DO 10 i=1,m sw=0 IF(leap.eq.0) GOTO 30 ycount=ycount+1 IF(ycount.ne.leap) GOTO 30 sw=1 ycount=0 leap=4 ! from now on every 4th year will be a leap year c####################################################################### c begin the simulation process c####################################################################### 30 continue pretmin=-100. pretmax=200. notdone=0 DO 20 j=1,idyr 40 u=ran3(iseed) IF(nw.eq.0) GOTO 60 c####################################################################### c if proceding day is dry, use a new generated random number c and do following: c####################################################################### IF(u.lt.p10(j)) THEN GOTO 90 ! u < p10(n), day n is dry ELSE GOTO 70 ! u >= p10(n), day n is wet END IF c####################################################################### c if proceding day is wet, use a new generated random number c and do following: c####################################################################### 60 IF(u.lt.p00(j)) GOTO 90 ! u < p00(n) day n is dry c####################################################################### c for wet day the mixed exponential distribution is: c####################################################################### 70 nw=1 u=ran3(iseed) ! new random number IF(u.lt.albar) GOTO 80 ! u < alpha(n): see orange book Eq.(18),Page 7 u=ran3(iseed) p=-delta(j)*alog(u)+0.01 GOTO 100 ! u >= alpha(n): see orange book Eq.(19),Page 7 80 u=ran3(iseed) p=-beta(j)*alog(u)+0.01 GOTO 100 c####################################################################### c for day is dry c####################################################################### 90 p=0. nw=0 GOTO 100 100 CONTINUE DO 120 mo=1,11 IF(j.ge.mday(mo).and.j.lt.mday(mo+1)) ji=mo 120 CONTINUE IF(j.ge.mday(12)) ji=12 mmd=j-mday(ji)+1 IF(notdone.eq.1) THEN mmd=mmd+1 notdone=0 END IF ijk=i if(ji.ge.11) ijk=i+1 c ccc IF(char1.eq.'Y'.or.char1.eq.'y') THEN ccc IF(char2.eq.'Y'.or.char2.eq.'y') THEN c CALL TMAX_TMIN_123(j,ns,iseed,nw,ji) ! simulate tmax,tmin,rad dew=dew-100. wind=abs(wind) c ccc ELSE ccc CALL TMAX_TMIN_RAD(j,ns,iseed,nw) ! simulate tmax,tmin,rad ccc END IF c tmax=tmax-100. tmin=tmin-100. tmin=min(tmin,pretmax) tmax=max(tmax,pretmin) pretmin=tmin pretmax=tmax c ccc END IF c 110 CONTINUE nmo = ji+2 IF(ji.gt.10) nmo = ji-10 ! nmo -- # of month c ccc IF(char2.eq.'Y'.or.char2.eq.'y') THEN c write(21,'(i3,i3,a5,i3,f6.2,5f11.2)')ijk,nmo,mname(ji), !mmd,p,tmax,tmin,dew,wind,rad if(ijk.eq.1.or.ijk.eq.(m+1))goto 159 if(mname(ji).eq.' Feb'.and.mmd.eq.29) goto 159 write(22,'(i3,i3,i3,f6.2,5f11.2)')ijk-1,nmo, !mmd,p,tmax,tmin,dew,wind,rad 159 continue c IF(j.lt.idyr) GOTO 20 IF(sw.eq.0) GOTO 20 sw=0 notdone=1 GOTO 40 20 CONTINUE 10 CONTINUE close(22) RETURN END c ################################################################## c ################################################################## c ###### ###### c ###### SUBROUTINE DISAGGR ###### c ###### ###### c ###### (c) 1994 ###### c ###### ###### c ################################################################## c ################################################################## c ################################################################## c ################################################################## c ###### ###### c ###### SUBROUTINE NCAR ###### c ###### ###### c ###### (c) 1994 ###### c ###### ###### c ################################################################## c ################################################################## c ################################################################## c ################################################################## c ###### ###### c ###### SUBROUTINE FILTER ###### c ###### ###### c ###### (c) 1995 ###### c ###### ###### c ################################################################## c ################################################################## c ################################################################## c ################################################################## c ###### ###### c ###### SUBROUTINE READ_ESTIMATION ###### c ###### ###### c ###### (c) 1994 ###### c ###### ###### c ################################################################## c ################################################################## SUBROUTINE READ_ESTIMATION(ii,stname,lat) c c####################################################################### c c PURPOSE: c c Read temperture and radiation parameters from file c "temp.dat". c Temperture and radiation procedure is from Richardson, C.W. and c D.A. Wright 'WGEN: A model for generating daily weather varibles' c USDA-ARS-8, 1984: 83pp. c c####################################################################### c c AUTHOR: C.L.Hanson,K.A.Cumming,D.A.Woolhiser and C.W.Richardson c c Modify history: c Yunyun Lu (November,1994) c transfer program from QBASIC to FORTRAN 77 c c Y. Lu (Dec. 1994) c added the full documentation c c####################################################################### c c INPUT: c c instna station's name. c w temperatures in F, radiation in Langleys c 1. Annual mean of tmax for dry days; c 2. Amplitude of tmax for wet or dry days; c 3. Annual mean of the coefficient of variation of tmax c for wet or dry days; c 4. Amplitude of the coefficient of variation of tmax c for wet or dry days; c 5. Annual mean of tmax for wet days; c 6. Annual mean of tmax for wet or dry days; c 7. Amplitude of tmin for wet or dry days; c 8. Annual mean of the coefficient of variation of tmin c for wet or dry days; c 9. Amplitude of the coefficient of variation of tmin c for wet or dry days; c 10. Annual mean of radiation for dry days; c 11. Amplitude of radiation for dry days; c 12. Annual mean of radiation for wet days; c 13. Amplitude of radiation for wet days; c 14. Annual mean of the coefficient of variation of c radiation for dry days; c 15. Amplitude of the coefficient of variation of c radiation for dry days; c 16. Annual mean of the coefficient of variation of c radiation for dry days; c 17. Amplitude of the coefficient of variation of c radiation for dry days; c c OUTPUT: c c tmad mean of tmax -- dry days c txs standard deviation of tmax -- dry days c txm1 mean of tmax -- wet days c txs1 standard deviation of tmax -- wet days c tnm mean of tmin -- wet and dry days c tns standard deviation of tmin -- wet and dry days c rm0 mean of radiation -- dry days c rs0 standard deviation of radiation -- dry days c rm1 mean of radiation -- wet days c rs1 standard deviation of radiation -- wet days c rm max. radiation along a given latitude c chipre zero array c c WORKING ARRAY: c c sw Averaged w by using stations' info within 100 miles c c####################################################################### c c Variable Declarations. c c####################################################################### c parameter(m=360,n=4,nh=13,idyr=365) common /const1/pi common /sited2/delta(idyr),p00(idyr),p10(idyr),beta(idyr),albar common /tem11/tmad(idyr),txs(idyr),txm1(idyr),txs1(idyr), !dmad(idyr),dxs(idyr),dxm1(idyr),dxs1(idyr), !wmad(idyr),wxs(idyr),wxm1(idyr),wxs1(idyr), !tnm(idyr),tns(idyr),tnm1(idyr),tns1(idyr), !rm0(idyr),rs0(idyr),rm1(idyr),rs1(idyr), !rm(idyr),chipre(6) real lat !latitude of station real w(nh) !input parameters *.temp real sw(17) !average of input parameters character*20 stname(ii),instna !station's name. character*25 site,temp,matrix,gem6out,filename common/filname/site,temp,matrix,gem6out,filename C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C C Beginning of executable code... C C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DO 10 k=1,nh sw(k)=0. 10 CONTINUE np=0 open(2,file='usclimate.temp',status='old') 15 read(2,'(A20)',END=30) instna ! read station name read(2,*) (w(i1),i1=1,nh) ! read parameter c####################################################################### c find the stations within 100 miles (ii>1) or closest station (ii=1) c####################################################################### DO 40 i=1,ii IF(instna.eq.stname(i)) THEN np=np+1 DO 20 i1=1,nh sw(i1)=sw(i1)+w(i1) 20 CONTINUE END IF 40 CONTINUE GOTO 15 30 close(2) DO 50 i=1,nh IF (np.ne.0) sw(i)=sw(i)/(float(np)) 50 CONTINUE c####################################################################### c calculate the max. radiation along lat. c####################################################################### DO 60 i=1,idyr sd=0.4102*sind((360/365.2425)*(i-21.25)) !corrected for solar day ch=-tand(lat)*tan(sd) IF (ch.gt.1.) ch=1. IF (ch.lt.-1.) ch=-1. h=acos(ch) dd=1.+0.0335*sind((360/365.2425)*(i+88.2)) !corrected for solar day c rm(i)=889.2305*dd*(h*sind(lat)*sin(sd)+cos(sd)*sin(h)) c ! ^ (cosd(lat)*cos(sd)... rm(i)=889.2305*dd*(h*sind(lat)*sin(sd)+(cosd(lat)*cos(sd)*sin(h))) c c rm(i)=rm(i)*0.9 ! 0.9 in Richardson paper rm(i)=rm(i)*0.8 60 CONTINUE c####################################################################### c add 100 degrees to temperature related parameter c####################################################################### sw(1)=sw(1)+100. sw(5)=sw(5)+100. sw(6)=sw(6)+100. c####################################################################### c change solar radiation CV parameters if necessary c####################################################################### sw(14)=.224 sw(15)=-.084 sw(16)=.488 sw(17)=-.135 c####################################################################### c calculate daily parameter values for tmax, tmin and radiation c####################################################################### DO 70 i=1,idyr c ccc dt=cosd(float(i-141)) ccc dr=cosd(float(i-113)) c !!! cos((2pi/365.2425)*float(i-141)) c ! cosd((360/365.2425)*float(i-141)) phase angle dt=cosd((360/365.2425)*float(i-141)) !correction for solar day dr=cosd((360/365.2425)*float(i-113)) !!! why 141 and 113 days? c tmad(i)=sw(1)+sw(2)*dt ! mean of tmax -- dry days xcr=sw(3)+sw(4)*dt txs(i)=tmad(i)*xcr ! std dev tmax -- dry days c txm1(i)=tmad(i)-sw(1)+sw(5) ! mean of tmax-wet days txs1(i)=txm1(i)*xcr ! std dev of tmax -- wet days c tnm(i)=sw(6)+sw(7)*dt ! mean of tmin--wet & dry days xcr=sw(8)+sw(9)*dt IF (xcr.lt.0.) xcr=0.06 ! WHY 0.06 ??? tns(i)=tnm(i)*xcr !std dev of tmin-wet & dry days c rm0(i)=sw(10)+sw(11)*dr ! mean radiation -- dry days xcr=sw(14)+sw(15)*dr rs0(i)=rm0(i)*xcr ! std dev rad -- dry days c rm1(i)=sw(12)+sw(13)*dr ! mean radiation -- wet days xcr=sw(16)+sw(17)*dr rs1(i)=rm1(i)*xcr !std dev rad -- wet days c 70 CONTINUE c####################################################################### c subtract 100 degrees to temperature related parameter c####################################################################### sw(1)=sw(1)-100. sw(5)=sw(5)-100. sw(6)=sw(6)-100. c####################################################################### c for the first run chi is zero. see orange book Eq. (13), Page 5 c####################################################################### DO 80 i=1,3 chipre(i)=0. 80 CONTINUE RETURN END c ################################################################## c ################################################################## c ###### ###### c ###### SUBROUTINE TMAX_TMIN_RAD ###### c ###### ###### c ###### (c) 1994 ###### c ###### ###### c ################################################################## c ################################################################## SUBROUTINE TMAX_TMIN_RAD(j,ns,iseed,nw) c c####################################################################### c c PURPOSE: c c This program simulates max and min daily temp and radiation c c####################################################################### c c AUTHOR: C.L.Hanson,K.A.Cumming,D.A.Woolhiser and C.W.Richardson c c Modify history: c Yunyun Lu (November,1994) c transfer program from QBASIC to FORTRAN 77 c c Y. Lu (Dec. 1994) c added the full documentation c c####################################################################### c c INPUT: c c tmad mean of tmax -- dry days c txs standard deviation of tmax -- dry days c txm1 mean of tmax -- wet days c txs1 standard deviation of tmax -- wet days c tnm mean of tmin -- wet and dry days c tns standard deviation of tmin -- wet and dry days c rm0 mean of radiation -- dry days c rs0 standard deviation of radiation -- dry days c rm1 mean of radiation -- wet days c rs1 standard deviation of radiation -- wet days c rm max. radiation along a given latitude c chipre last step of chi(n), =chi(n-1) c c OUTPUT: c c p Simulated M years of daily rainfall c tmax Simulated M years of daily max. temperature c tmin Simulated M years of daily min. temperature c rad Simulated M years of daily total radiation c c####################################################################### c c Variable Declarations. c c####################################################################### c parameter(idyr=365) common /const1/pi common /tem11/tmad(idyr),txs(idyr),txm1(idyr),txs1(idyr), !dmad(idyr),dxs(idyr),dxm1(idyr),dxs1(idyr), !wmad(idyr),wxs(idyr),wxm1(idyr),wxs1(idyr), !tnm(idyr),tns(idyr),tnm1(idyr),tns1(idyr), !rm0(idyr),rs0(idyr),rm1(idyr),rs1(idyr), !rm(idyr),chipre(6) common /temp2/p,tmax,tmin,dew,wind,rad common /temp3/fran(12) real a(3,3),b(3,3) data a/.567,.253,-.006,.086,.504,-.039,-.002,-.05,.244/ data b/.781,.328,.238,.0,.637,-.341,.0,.0,.873/ real rd(3),rrd(3),epsilon(3),chi(3) character*25 site,temp,matrix,gem6out,filename common/filname/site,temp,matrix,gem6out,filename C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C C Beginning of executable code... C C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% rmm=rm0(j) rs=rs0(j) txxm=tmad(j) txxs=txs(j) IF (nw.eq.1) THEN rmm=rm1(j) rs=rs1(j) txxm=txm1(j) txxs=txs1(j) END IF IF(ns.le.1) THEN DO k=1,5,2 u=ran3(iseed) u1=ran3(iseed) fran(k)=sqrt(-2*alog(u))*cos(2*pi*u1) fran(1+k)=sqrt(-2*alog(u))*sin(2*pi*u1) END DO ns=3 ELSE ns=0 END IF c####################################################################### c for first day using last three number in the fran() array, THEN c for next day using the first three number in fran() array. c Iterately do it back and forewards. c integer ns control which three number will be used. c####################################################################### DO k=1,3 epsilon(k)=fran(k+ns) rd(k)=0. rrd(k)=0. END DO DO k=1,3 DO l=1,3 rd(k)=rd(k)+b(k,l)*epsilon(l) ! see orange book Eq.(13) page 5 rrd(k)=rrd(k)+a(k,l)*chipre(l) END DO END DO DO l=1,3 chi(l)=rd(l)+rrd(l) chipre(l)=chi(l) ! at the END of calculation ! store chi as next run's chipre. END DO tmax=chi(1)*txxs+txxm tmin=chi(2)*tns(j)+tnm(j) c####################################################################### c Adjust for tmax, tmin, radiation c####################################################################### IF(tmin.gt.tmax) THEN ! tmax should never be smaller than tmin tmaxt=tmax tmax=tmin tmin=tmaxt END IF rad=chi(3)*rs+rmm ! rad will not smaller than 0.1*max.radation rmin=.1*rm(j) ! and not larger than max. radiation rad=max(rad,rmin) rad=min(rad,rm(j)) RETURN END c ################################################################## c ################################################################## c ###### ###### c ###### SUBROUTINE READ_123IMATION ###### c ###### ###### c ###### (c) 1994 ###### c ###### ###### c ################################################################## c ################################################################## SUBROUTINE READ_123IMATION(lat) c c####################################################################### c c PURPOSE: c c Read temperture and radiation parameters from file c "temp.dat". c Temperture and radiation procedure is from Richardson, C.W. and c D.A. Wright 'WGEN: A model for generating daily weather varibles' c USDA-ARS-8, 1984: 83pp. c c####################################################################### c c AUTHOR: C.L.Hanson,K.A.Cumming,D.A.Woolhiser and C.W.Richardson c c Modify history: c Yunyun Lu (November,1994) c transfer program from QBASIC to FORTRAN 77 c c Y. Lu (Dec. 1994) c added the full documentation c c####################################################################### c c INPUT: c c instna station's name. c w temperatures in F, radiation in Langleys c 1. Annual mean of tmax for dry days; c 2. Amplitude of tmax for wet or dry days; c 3. Annual mean of the coefficient of variation of tmax c for wet or dry days; c 4. Amplitude of the coefficient of variation of tmax c for wet or dry days; c 5. Annual mean of tmax for wet days; c 6. Annual mean of tmax for wet or dry days; c 7. Amplitude of tmin for wet or dry days; c 8. Annual mean of the coefficient of variation of tmin c for wet or dry days; c 9. Amplitude of the coefficient of variation of tmin c for wet or dry days; c 10. Annual mean of radiation for dry days; c 11. Amplitude of radiation for dry days; c 12. Annual mean of radiation for wet days; c 13. Amplitude of radiation for wet days; c 14. Annual mean of the coefficient of variation of !?? see 16 c radiation for dry days; c 15. Amplitude of the coefficient of variation of !?? see 17 c radiation for dry days; c 16. Annual mean of the coefficient of variation of c radiation for dry days; c 17. Amplitude of the coefficient of variation of c radiation for dry days; c c OUTPUT: c c tmad mean of tmax -- dry days c txs standard deviation of tmax -- dry days c txm1 mean of tmax -- wet days c txs1 standard deviation of tmax -- wet days c tnm mean of tmin -- wet and dry days c tns standard deviation of tmin -- wet and dry days c rm0 mean of radiation -- dry days c rs0 standard deviation of radiation -- dry days c rm1 mean of radiation -- wet days c rs1 standard deviation of radiation -- wet days c rm max. radiation along a given latitude c chipre zero array c c WORKING ARRAY: c c sw Averaged w by using stations' info within 100 miles c c####################################################################### c c Variable Declarations. c c####################################################################### c parameter(m=360,n=4,nh=40,idyr=365) common /const1/pi common /sited2/delta(idyr),p00(idyr),p10(idyr),beta(idyr),albar common /tem11/tmad(idyr),txs(idyr),txm1(idyr),txs1(idyr), !dmad(idyr),dxs(idyr),dxm1(idyr),dxs1(idyr), !wmad(idyr),wxs(idyr),wxm1(idyr),wxs1(idyr), !tnm(idyr),tns(idyr),tnm1(idyr),tns1(idyr), !rm0(idyr),rs0(idyr),rm1(idyr),rs1(idyr), !rm(idyr),chipre(6) common /temp4/a(6,6,12),b(6,6,12) real lat real sw(nh) !average of input parameters character*25 site,temp,matrix,gem6out,filename common/filname/site,temp,matrix,gem6out,filename C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C C Beginning of executable code... C C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DO 10 k=1,nh sw(k)=0. 10 CONTINUE np=0 ccc open(2,file='boise61a.temp',status='old') open(2,file=temp,status='old') read(2,*) (sw(i1),i1=1,nh) ! read annual means parameters c####################################################################### c find the stations within 100 miles (ii>1) or closest station (ii=1) c####################################################################### 30 close(2) c####################################################################### c calculate the max. radiation along lat. c####################################################################### DO 60 i=1,idyr sd=0.4102*sind((360/365.2425)*(i-21.25)) !corrects to solar day ch=-tand(lat)*tan(sd) ch=min(ch,1.) ch=max(ch,-1.) h=acos(ch) dd=1.+0.0335*sind((360/365.2425)*(i+88.2)) ! ccc rm(i)=889.2305*dd*(h*sind(lat)*sin(sd)+cos(sd)*sin(h)) c ! ^ (cosd(lat)*cos(sd)... rm(i)=889.2305*dd*(h*sind(lat)*sin(sd)+(cosd(lat)*cos(sd)*sin(h))) c !REF WGEN & WGENPAR c rm(i)=rm(i)*0.9 rm(i)=rm(i)*0.8 !Ref WGEN & WGENPAR c 60 CONTINUE c c####################################################################### c add 100 degrees to temperature related parameter c####################################################################### c do io=1,21,4 sw(io)=sw(io)+100. end do c####################################################################### c change solar radiation CV parameters if necessary c####################################################################### c fixed in US Climate ARS-114 July 1994, c Hanson, Cumming, Woolhiser, Richardson c sw(35)=.224 c sw(36)=-.084 c sw(39)=.488 c sw(40)=-.135 c####################################################################### c calculate daily parameter values for tmax, tmin and radiation c####################################################################### DO 70 i=1,idyr c !360/365.2425 deg / solar day dt=cosd((360/365.2425)*float(i-141)) ! ?? 141,35,30,113 dw=cosd((360/365.2425)*float(i-35-30)) dr=cosd((360/365.2425)*float(i-113)) c tmad(i)=sw(1)+sw(2)*dt ! mean of tmax -- dry days xcr=sw(3)+sw(4)*dt txs(i)=tmad(i)*xcr ! standard deviation of tmax -- dry days c txm1(i)=sw(5)+sw(6)*dt ! mean of tmax-wet days xcr=sw(7)+sw(8)*dt txs1(i)=txm1(i)*xcr ! standard deviation of tmax -- wet days c tnm(i)=sw(9)+sw(10)*dt ! mean of tmin -- dry days xcr=sw(11)+sw(12)*dt c tns(i)=tnm(i)*xcr ! standard deviation of tmin -- dry days tnm1(i)=sw(13)+sw(14)*dt ! mean of tmin -- wet days xcr=sw(15)+sw(16)*dt tns1(i)=tnm1(i)*xcr ! standard deviation of tmin -- wet days c dmad(i)=sw(17)+sw(18)*dt ! mean of dew -- dry days xcr=sw(19)+sw(20)*dt dxs(i)=dmad(i)*xcr ! standard deviation of dew -- dry days c dxm1(i)=sw(21)+sw(22)*dt ! mean of dew-wet days xcr=sw(23)+sw(24)*dt dxs1(i)=dxm1(i)*xcr ! standard deviation of dew -- wet days c rm0(i)=sw(25)+sw(26)*dr ! mean radiation -- dry days xcr=sw(27)+sw(28)*dr rs0(i)=rm0(i)*xcr ! standard deviation of radiation -- dry days c rm1(i)=sw(29)+sw(30)*dr ! mean radiation -- wet days xcr=sw(31)+sw(32)*dr rs1(i)=rm1(i)*xcr !standard deviation of radiation -- wet days c wmad(i)=sw(33)+sw(34)*dw ! mean of wind -- dry days xcr=sw(35)+sw(36)*dw wxs(i)=wmad(i)*xcr ! standard deviation of wind -- dry days c wxm1(i)=sw(37)+sw(38)*dw ! mean of wind-wet days xcr=sw(39)+sw(40)*dw wxs1(i)=wxm1(i)*xcr ! standard deviation of wind -- wet days c 70 CONTINUE c####################################################################### c subtract 100 degrees to temperature related parameter c####################################################################### do io=1,21,4 sw(io)=sw(io)-100. end do c####################################################################### c for the first run chi is zero. see orange book Eq. (13), Page 5 c####################################################################### DO 80 i=1,6 chipre(i)=0. 80 CONTINUE ccc open(unit=2,file='boise61a.max',status='old') !was boi_m_6.mat open(unit=2,file=matrix,status='old') !was boi_m_6.mat read(2,*) !win98 did not like .mat do jj=1,12 read(2,*) read(2,*) read(2,'(6f7.3)') ((a(i,j,jj),j=1,6),i=1,6) end do read(2,*) do jj=1,12 read(2,*) read(2,*) read(2,'(6f7.3)') ((b(i,j,jj),j=1,6),i=1,6) end do close(2) RETURN END c ################################################################## c ################################################################## c ###### ###### c ###### SUBROUTINE TMAX_TMIN_123 ###### c ###### ###### c ###### (c) 1994 ###### c ###### ###### c ################################################################## c ################################################################## SUBROUTINE TMAX_TMIN_123(j,ns,iseed,nw,mo) c c####################################################################### c c PURPOSE: c c This program simulates max and min daily temp and radiation c c####################################################################### c c AUTHOR: C.L.Hanson,K.A.Cumming,D.A.Woolhiser and C.W.Richardson c c Modify history: c Yunyun Lu (November,1994) c transfer program from QBASIC to FORTRAN 77 c c Y. Lu (Dec. 1994) c added the full documentation c c####################################################################### c c INPUT: c c tmad mean of tmax -- dry days c txs standard deviation of tmax -- dry days c txm1 mean of tmax -- wet days c txs1 standard deviation of tmax -- wet days c tnm mean of tmin -- wet and dry days c tns standard deviation of tmin -- wet and dry days c rm0 mean of radiation -- dry days c rs0 standard deviation of radiation -- dry days c rm1 mean of radiation -- wet days c rs1 standard deviation of radiation -- wet days c rm max. radiation along a given latitude c chipre last step of chi(n), =chi(n-1) c c OUTPUT: c c p Simulated M years of daily rainfall c tmax Simulated M years of daily max. temperature c tmin Simulated M years of daily min. temperature c rad Simulated M years of daily total radiation c c####################################################################### c c Variable Declarations. c c####################################################################### c parameter(idyr=365,nn=6) common /const1/pi common /tem11/tmad(idyr),txs(idyr),txm1(idyr),txs1(idyr), !dmad(idyr),dxs(idyr),dxm1(idyr),dxs1(idyr), !wmad(idyr),wxs(idyr),wxm1(idyr),wxs1(idyr), !tnm(idyr),tns(idyr),tnm1(idyr),tns1(idyr), !rm0(idyr),rs0(idyr),rm1(idyr),rs1(idyr), !rm(idyr),chipre(nn) common /temp2/p,tmax,tmin,dew,wind,rad common /temp3/fran(12) common /temp4/a(nn,nn,12),b(nn,nn,12) real rd(nn),rrd(nn),epsilon(nn),chi(nn) character*25 site,temp,matrix,gem6out,filename common/filname/site,temp,matrix,gem6out,filename c C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C C Beginning of executable code... C C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% rmm=rm0(j) rs=rs0(j) txxm=tmad(j) txxs=txs(j) tnnm=tnm(j) tnns=tns(j) dxxm=dmad(j) dxxs=dxs(j) wxxm=wmad(j) wxxs=wxs(j) IF (nw.eq.1) THEN rmm=rm1(j) rs=rs1(j) txxm=txm1(j) txxs=txs1(j) tnnm=tnm1(j) tnns=tns1(j) dxxm=dxm1(j) dxxs=dxs1(j) wxxm=wxm1(j) wxxs=wxs1(j) END IF IF(ns.le.1) THEN DO k=1,12,2 u=ran3(iseed) u1=ran3(iseed) fran(k)=sqrt(-2*alog(u))*cos(2*pi*u1) fran(1+k)=sqrt(-2*alog(u))*sin(2*pi*u1) END DO ns=nn ELSE ns=0 END IF c####################################################################### c for first day using last three number in the fran() array, THEN c for next day using the first three number in fran() array. c Iterately do it back and forewards. c integer ns control which three number will be used. c####################################################################### DO k=1,nn epsilon(k)=fran(k+ns) rd(k)=0. rrd(k)=0. END DO DO k=1,nn DO l=1,nn rd(k)=rd(k)+b(k,l,mo)*epsilon(l) ! see orange book Eq.(13) page 5 rrd(k)=rrd(k)+a(k,l,mo)*chipre(l) END DO END DO DO l=1,nn chi(l)=rd(l)+rrd(l) chipre(l)=chi(l) ! at the END of calculation ! store chi as next run's chipre. END DO tmax=chi(2)*txxs+txxm tmin=chi(3)*tnns+tnnm dew=chi(5)*dxxs+dxxm wind=chi(6)*wxxs+wxxm c print*,'wind',chi(5),wxxs,chi(5)*wxxs,wxxm,wind c####################################################################### c Adjust for tmax, tmin, radiation c####################################################################### IF(tmin.gt.tmax) THEN ! tmax should never be smaller than tmin tmaxt=tmax tmax=tmin tmin=tmaxt END IF rad=chi(4)*rs+rmm ! rad not smaller than 0.1*max.radation rmin=.1*rm(j) ! and not larger than max. radiation rad=max(rad,rmin) rad=min(rad,rm(j)) RETURN END c ################################################################## c ################################################################## c ###### ###### c ###### SUBROUTINE TOTAL_RAINFALL ###### c ###### ###### c ###### (c) 1994 ###### c ###### ###### c ################################################################## c ################################################################## SUBROUTINE TOTAL_RAINFALL c c####################################################################### c c PURPOSE: c c Calculate the cdf for total rainfall in n days. c Using equation from ToDOrovic, P. and D.A. Woolhiser. 1975 c A stochastic model of n-day precipitation. J.Applied Meteor. c 14(1): 17-24. xlam=1/(mean daily rainfall) c c####################################################################### c c AUTHOR: C.L.Hanson,K.A.Cumming,D.A.Woolhiser and C.W.Richardson c c Modify history: c Yunyun Lu (November,1994) c transfer program from QBASIC to FORTRAN 77 c c Y. Lu (Dec. 1994) c added the full documentation c c####################################################################### c c c INPUT: c c p00 Transition probability of a dry day following a dry day c p10 Transition probability of a dry day following a wet day c beta The mean of the smaller exponential distributions c delta The mean of the larger exponential distributions c albar The weighting function c c####################################################################### c c Variable Declarations. c c####################################################################### c parameter (idyr=365) parameter (n1=4,nh=13,myr=12) integer mday(myr) integer xd,xm common /const1/pi common /sited2/delta(idyr),p00(idyr),p10(idyr),beta(idyr),albar real xa(40),fx(40),psii(40),psi0(40) real exp character*20 char1 data mday/307,338,1,32,62,93,123,154,185,215,246,276/ ! start from March 1 !!! C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C C Beginning of executable code... C C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c print*,'Enter the month: (1,2,...,12)' read(5,'(I4)') xm print*,'Enter the day: (1,2,...,31)' read(5,'(I4)') xd print*, ! 'Enter the length,n,of the period in days (1,2,...,31)' read(5,'(I4)') n j=mday(xm)+xd-1+n/2 IF (j.gt.idyr) j=j-idyr c####################################################################### c calculate counting function density functions c####################################################################### q0=1-p00(j) qi=1-p10(j) psii(1)=(1-q0)**(n-1)*(1-qi) psi0(1)=(1-q0)**n nu=0 nt=n+1 DO 100 i=2,n nu=nu+1 nci=aint(n+0.5-abs(2*nu-n+0.5)+0.01) nc0=aint(n+0.5-abs(2*nu-n-0.5)+0.01) nc=0 a=0 bc=1 nsw=1 sumi=0 sum0=0 termi=(1-qi)/(1-q0) term0=q0/qi 10 sumi=sumi+termi sum0=sum0+term0 nc=nc+1 IF(nc.ne.nc0) GOTO 30 psi0(i)=qi**nu*(1-q0)**(n-nu)*sum0 IF(nc0.gt.nci) GOTO 100 IF(nc.eq.nci) GOTO 40 20 IF (nsw.eq.2) GOTO 50 nsw=2 a=a+1 termi=termi*(nu-a+1)/a*q0/qi term0=term0*(n-nu-a+1)/a*(1-qi)/(1-q0) GOTO 10 30 IF(nc.ne.nci) GOTO 20 40 psii(i)=qi**nu*(1-q0)**(n-nu)*sumi IF(nc0.gt.nci) GOTO 20 GOTO 100 50 nsw=1 bc=bc+1 termi=termi*(n-nu-bc+1)/(bc-1)*(1-qi)/(1-q0) term0=term0*(nu-bc+1)/(bc-1)*q0/qi GOTO 10 100 CONTINUE psii(nt)=qi**n psi0(nt)=qi**n*q0/qi c####################################################################### c in this section we computer the CDF for n days c set exponential parameter to 1/mean c####################################################################### xlam=1/(albar*beta(j)+(1-albar)*delta(j)) nc=0 xa(1)=0. k=j-n/2 IF(k.lt.1) k=idyr ra=(1-p00(k))/(1+p10(k)-p00(k)) sbar=ra*n/xlam dx=sbar/4 IF(dx.le.0.05) dx=0.05 print*, :'Was it precipitating on the day before', :xm,'-',xd,'? (Y,N,E(don"t know))' read(5,'(a1)') char1 IF(char1.eq.'Y'.or.char1.eq.'y') ra=1. IF(char1.eq.'N'.or.char1.eq.'n') ra=0. IF(char1.eq.'E'.or.char1.eq.'e') THEN write(6,'(/ 22(/5x,a)//)') :'Enter the probability of precipitation :on the day before in the form 0.xx' read(5,'(f10.5)')rb IF(rb.ne.0.0) ra=rb END IF fx(1)=(1-q0)-(ra*(qi-q0)) fx(1)=fx(1)*(1-q0)**(n-1) sav=fx(1) DO 110 i=2,40 xa(i)=xa(i-1)+dx sum=0.0 coef=xlam nc=nc+1 DO 120 nu=1,n co=xa(i) suma=0.0 terma=xlam*co**nu/nu DO 130 j=1,nu terma=terma*(nu-j+1)/(xlam*co) suma=suma-terma 130 CONTINUE xc=terma/xlam+suma/(exp(xlam*co)*xlam) term=(ra*psii(nu+1)+(1-ra)*psi0(nu+1))*coef*xc sum=sum+term coef=coef*xlam/nu fx(i)=sav+sum 120 CONTINUE IF (fx(i).lt.0.99) GOTO 110 in=i GOTO 140 110 CONTINUE in=40 140 print*, !'Cumulative distribution function for',n,'days.' print*,'Begining',xm,'-',xd IF(ra.eq.0.0.or.ra.eq.1.0) THEN IF(ra.eq.1) char1='WET' IF(ra.eq.0) char1='DRY' print*,'Day before ',n,' - day period was ',char1 END IF IF(ra.gt.0.0.or.ra.lt.1.0) THEN print*,' x (inches) f(x)' DO 150 i=1,in write(6,'(2f15.8)') xa(i),fx(i) 150 CONTINUE END IF RETURN END c ################################################################## c ################################################################## c ###### ###### c ###### SUBROUTINE DISTANCE_SORT ###### c ###### ###### c ###### (c) 1994 ###### c ###### ###### c ################################################################## c ################################################################## SUBROUTINE DISTANCE_SORT(stname,z,zz,ii,r) c c####################################################################### c c PURPOSE: c c Put data in order: from max. to min. or from min. to max. c c####################################################################### c c AUTHOR: Yunyun Lu c c####################################################################### c c Variable Declarations. c c####################################################################### c PARAMETER (n=4,nh=13) real r(ii),z(ii,n,nh),zz(ii),a character*20 stname(ii),b c####################################################################### c c C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C C Beginning of executable code... C C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DO 10 I=1,II-1 ISMALL=I DO 20 J=I+1,II IF(R(J).LT.R(ISMALL)) ISMALL=J 20 CONTINUE DO 30 I1=1,N DO 30 I2=1,NH A=Z(I,I1,I2) Z(I,I1,I2)=Z(ISMALL,I1,I2) Z(ISMALL,I1,I2)=A 30 CONTINUE B=STNAME(I) STNAME(I)=STNAME(ISMALL) STNAME(ISMALL)=B A=ZZ(I) ZZ(I)=ZZ(ISMALL) ZZ(ISMALL)=A 10 CONTINUE RETURN END c c ################################################################## c ################################################################## c ###### ###### c ###### FUNCTION RAN3 ###### c ###### ###### c ###### Copyright (c) 1995 ###### c ###### ###### c ################################################################## c ################################################################## c FUNCTION RAN3(iseed) c c####################################################################### c c PURPOSE: c c Generates a random number between 0 and 1 by feeding c a negative integer iseed. c c Reference: "Seminumerical Algorithms" by Donald Knuth c c####################################################################### c c INPUT : c c iseed an arbitrary negative integer as a seed for a c sequence of random numbers c c OUTPUT: c c ran3 A random number between 0 and 1. c c####################################################################### c c c####################################################################### c c Variable Declarations: c c####################################################################### c implicit none ! Force explicit declarations integer iseed ! The seed for random number generation real ran3 ! The function to generate random number. c c####################################################################### c c Miscellaneous local variables: c c####################################################################### c integer mbig,mseed,mz,k,ii,inext,inextp,i,iff,mk,mj real fac parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1.e-9) integer ma(55) save iff, inext, inextp, ma data iff /0/ c C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C C Beginning of executable code... C C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c c c###################################################################### c c Initialize the sequence of random numbers between 0 and 1, c using iseed. c c###################################################################### c IF (iseed.lt.0.or.iff.eq.0) THEN iff=1 mj=mseed-iabs(iseed) mj=mod(mj,mbig) ma(55)=mj mk=1 DO 11 i=1,54 ii=mod(21*i,55) ma(ii)=mk mk=mj-mk IF (mk.lt.mz) mk=mk+mbig mj=ma(ii) 11 CONTINUE DO 13 k=1,4 DO 12 i=1,55 ma(i)=ma(i)-ma(1+mod(i+30,55)) IF (ma(i).lt.mz) ma(i)=ma(i)+mbig 12 CONTINUE 13 CONTINUE inext=0 inextp=31 iseed=1 ENDIF c c###################################################################### c c Start to generate a random number. c c###################################################################### c inext=inext+1 IF (inext.eq.56) inext=1 inextp=inextp+1 IF (inextp.eq.56) inextp=1 mj=ma(inext)-ma(inextp) IF (mj.lt.mz) mj=mj+mbig ma(inext)=mj ran3=mj*fac RETURN END