! This version has been updated to compile with the Lahey lf95 PC compiler ! It has not been tested under Unix or with the Watcom PC compiler since ! then (this would require a minimum of switching the compiler specific ! code to read in commandline args). It is possible that the random ! seed determined from the processor clock would need to be changed ! as well - there used to be compiler specific code there, but it isn't ! in this copy anymore. ! Wed May 17 01:11:34 CDT 2000 - LEW c Program Cligen. WEPP Water Erosion Project Durant, OK. c Please address inquiries to c c Arlin Nicks c USDA- ARS c PO Box 1430 c Durant, OK 74702 c Phone 405 924 5066 c c Version 2.1 October 1990 Durant, OK. c Multiple year generation, output file, and summary added c Dewpoint temperature added with wind speed and velocity c using mean standard deviation and sqew coefficient. c Change made to place parameter in state files (i.e. al Alabama c parms). All parameters rainfall, temperature, radiation, dewpoint c temperature, wind speed and direction are combined into one set c with 82 lines of data. c cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX c c Mon May 20 15:53:40 CDT 1991 - Rob Quinn c If a 'configure' file exists, read parameters from it c rather than being interactive. c No longer implemented - LEW c c Wed Jun 5 16:29:14 CDT 1991 - Rob Quinn c Read command line arguments. c c Tue May 16 23:16:30 CDT 2000 - LEW c Read command line arguments with Lahey lf95 PC compiler. c Code for reading command line args with Unix compiler c commented out using a "C" and labeled as such. c Code for reading command line args with Watcom PC compiler c commented out using a "*" and labeled as such. c c Code last updated: Tue May 16 23:04:36 CDT 2000 - LEW c cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX c common /bk1/ wvl(16,4,12),dir(12,17),wv,th,pi2,ang,rh(12), 1 tdp,calm(12) common /bk3/ alai(5),salb(5),dm(5),rd(5),uw,ep,es,rain,ws, 1 je,j,ida common /bk4/ nc(13),iyr,nt,mo,px,iopt,dtp(4),dmxi(4) common /bk7/ rst(12,3),t0(5),prw(2,12),obmx(12),obmn(12), 1 obsl(12),wft(12),subp(5),sno(5),snoev,k1(4), 1 k2(4),cvs(12),k3(4),k4(4),k5(4),k6(4),k7(4), 1 k8(4),k9(4),cvtx(12),cvtm(12),v1,v3,v5,v7,v9,v11,ra, 1 tmxg,tmng,rmx,yls,ylc,pit,st0,amp,avt,ts, 1 nsim,msim,l,tmx(5),tmn(5),tx(5),tb(5),to(5), 1 dmm(12,5),dm1(5),stdtx(12),stdtm(12),stdsl(12) c prw(1,12) monthly probability of wet day after dry day c prw(2,12) monthly probability of wet day after wet day common /bk5/ r(366),vmx(5),v(5),flu(6),bp(5),hc(5),hcr(5),fp(5), 1 cs(5),sk(5),smx(5),ev,sp,rl,q1,o,qd,af,pq,yp,aff, 1 amps,cn,sml,cfp(5),swm(5),sfc(5),dk(5) common /bk9/ hyd(250,5),xi(2,5),wi(12),tc(6),al(5),ab,ab1,rn1,r1, 1 pr,lu1,vo,nmb,tp,tr,ddt,tt common /cl1/ prcip(12,31),tgmx(12,31),tgmn(12,31),radg(12,31), 1 dur(12,31) common /sumr/sump(13),sumtx(12),sumtm(12),sumrd(12), 1 sumdr(12) c dimension sumpp(13),sumptx(12),sumptm(12),sumprd(12), 1 sumpdr(12) c dimension sm(100),smy(12),wind(366),windd(366),wgt(3),tymax(4) c ! Why is the next line commented out? Tue May 16 23:13:17 CDT 2000 - LEW c integer wlatd,wlatm,wlond,wlonm,elev c character*19 site(3) character*2 dr(54) character*1 yc c character*6 nstat !this commented out because it is an integer - LEW integer nstat !now declared as an integer - LEW character*41 stid character*41 stidd character*128 outfil !cligen db input file name(interactive mode) character*128 infile !cligen output file name character*25 twnam character*20 county cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX c Additional variables for non-interactive mode include 'f2kcli.inc' !For Fortran 2k commandline parsing code * include 'fsublib.fi' !Watcom integer function narg !Lahey character*128 cmdline !Lahey character*128 option !commandline option/arg str character*128 inputfile !cligen db input filename logical interactive, doheader !flags integer mxargs !maximum cmdline args integer argcnt !count number of cmdline args integer istate !cmdline arg for state no. integer ti(4) !variable for random seeds integer arg_len, arg_stat !Fortran 2k cmdline call arguments cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX c John Tatarko 07/16/93 c - compute an average annual air temperature real mairtmp, sumtmp c John Tatarko 02/02/98 c - compute monthly average air temperature (oC) real motemp(12) cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX data dr/'al','az','ar','ca','co','ct','de','fl','ga','id', 1 'il','in','ia','ks','ky','la','me','md','ma','mi', 1 'mn','ms','mo','mt','ne','nv','nh','nj','nm','ny', 1 'nc','nd','oh','ok','or','pa','ri','sc','sd','tn', 1 'tx','ut','vt','va','wa','wv','wi','wy','dc','ak', 1 'hi','pr','pi',' '/ data tymax/180.34,154.94,307.34,330.2/ c cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ! initial values for non-interactive code - LEW interactive=.true. doheader=.true. istate=-1 istatno= -1 ibyear=-1 numyr=-1 outfil='XXX' inputfile='XXX' cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX cX Old code for reading a parameter file - LEW cX If the file does not exist, continue as normal cX open(11,file='cligen.prm',err=07,status='old') cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ! LEW - Tue May 16 09:55:16 CDT 2000 ! determine if commandline arguments are specified ! if not , then use interactive mode C mxargs = nargs() !Unix * mxargs = iargc() !Watcom C mxargs = narg() !Lahey mxargs = command_argument_count() !Fortran 2k standard ! write(6,*) 'mxargs: ', mxargs if (mxargs .eq. 1) then !run in interactive mode goto 07 else !run in non-interactive mode interactive=.false. doheader=.false. endif cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX cX Old code for reading a parameter file - LEW cX01 read(11,'(1a)',end=05,err=05) option cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ! LEW - Tue May 16 09:55:16 CDT 2000 ! All commandline arguments are deciphered here ! both Unix and Watcom have a function to determine ! the number of arguments on the commandline. ! Lahey doesn't, so we must use "index() ! to determine the number of args ourselves first. !Changed for Lahey - don't know if other compilers need it to start at 0 !Wed May 17 01:01:17 CDT 2000 - LEW argcnt=1 01 argcnt=argcnt+1 !skip command name if(argcnt .gt. mxargs) goto 05 !no more args? !get arg/option C call getarg(argcnt,option,i) !Unix * i=igetarg(argcnt,option) !Watcom C call argopt(argcnt,option) !Lahey !Fortran 2k standard call get_command_argument(argcnt, option, arg_len, arg_stat) c State Number if (option(1:2).eq.'-S') then read(unit=option(3:),fmt='(1I2)') istate cXXXXXXXXXXXXXXXX above 14 lines modified by wagner and tatarko c 24/02/1998 XXXXXXXXXXXXXXXXXXXXXXXXXXXXX C write(6,*) 'State code:',istate c Station Number else if (option(1:2).eq.'-s') then read(option(3:),*) istatno C write(*,*) 'istatno:',istatno c Index Position else if(option(1:2).eq.'-x') then !not implemented continue c Random Seed value else if(option(1:2).eq.'-r') then if(option(3:).eq.'0') then !use default seed values continue else if(option(3:4).eq.'-1') then !use totally random seeds call gettim(ti(1),ti(2),ti(3),ti(4)) ti(3)=ti(4)+ti(3)*7+ti(2)*5+ti(1)*3 ti(3)=ti(3)-(ti(3)/199)*199 if (ti(3).lt.0) ti(3)=0-ti(3) c write(*,*) 'Advancing:',ti(3) do 2 j=1,ti(3) fx=randn(k1) fx=randn(k2) fx=randn(k3) fx=randn(k4) fx=randn(k5) fx=randn(k6) fx=randn(k7) fx=randn(k8) 2 fx=randn(k9) else !use specified seed value read(option(3:),*) i do 3 j=1,i fx=randn(k1) fx=randn(k2) fx=randn(k3) fx=randn(k4) fx=randn(k5) fx=randn(k6) fx=randn(k7) fx=randn(k8) 3 fx=randn(k9) endif c Beginning Year else if(option(1:2).eq.'-b') then read(option(3:),*) ibyear c write(*,*)'Beginning year:',ibyear c Duration in years else if(option(1:2).eq.'-y') then read(option(3:),*) numyr c write(*,*) 'Duration in years:',numyr c Hourly flag else if(option(1:2).eq.'-h') then continue !not implemented c Output file else if(option(1:2).eq.'-o') then outfil=trim(option(3:)) c write(*,*) 'Output file name:',outfil c Input File else if(option(1:2).eq.'-i') then inputfile=trim(option(3:)) c Threshold else if(option(1:2).eq.'-t') then !not implemented continue c Header flag else if(option(1:2).eq.'-H') then doheader=.true. c Verbose flag else if(option(1:2).eq.'-v'.or.option(1:2).eq.'-V') then write(*,*) 'Cligen - climate generator V2.1 Oct. 1990' write(*,*) 'Modified for non-interactive use' write(*,*) 'First modifications: Tuesday May 21 91' write(*,*) 'Last updated Wed Apr 4 17:24:48 2001' write(*,*) 'USDA-ARS Wind Erosion Research Unit' write(*,*) 'Throckmorton Hall, Kansas State University' c Cmdline syntax help flag else if(option(1:2).eq.'-?')then write(*,*) 'Cligen V2.1 - climate generator' write(*,*) 'Usage:' write(*,*) ' cligen -S -s' write(*,*) ' -i -o' write(*,*) ' -y -b' write(*,*) ' [-H] [-r] [-v] [-V] [-?]' write(*,*) 'Make sure that you use no spaces between each flag' write(*,*) 'and its parameter.' write(*,*) 'If no command line options are given, cligen' write(*,*) 'will interactively request the required info,' write(*,*) 'If any required comand line options are given,' write(*,*) 'then all other required command line options' write(*,*) 'must be specified as well.' write(*,*) 'Only -H, -V, -v, -r, and -? are optional.' stop else write(*,*) 'Unknown option:', option write(*,*) ' ignored' endif goto 01 cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX c05 close(unit=11,err=07) cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 05 goto 09 07 continue cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX call header write(*,*)' Press Enter to Continue' read(*,'(a1)')yc cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX c rjq 09 if (interactive.or.((istate.eq.(-1).and.inputfile.eq.'XXX'))) then go to 10 endif go to 51 cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 10 write(*,380) write(*,360) write(*,370) write(*,700) read(*,390,err=15)istate go to 16 15 write(*,*)' Error Entering State Code ' write(*,*) write(*,*)' Press Enter to Continue' read(*,'(a1)')yc go to 10 16 open(11,file='stations',status='old') rewind (11) iscnt=0 iyr=0 write(*,710) kknt=0 20 read(11,410,end=30)stid,nstate if(nstate.eq.istate) then kknt=kknt+1 iscnt=iscnt+1 backspace(11) read(11,420)stid,nst,nstat,lat,long,county if(kknt.gt.20) then kknt=0 write(*,*) 'Press "ENTER" to continue.' read(*,'(a1)') yc endif write(*,430)stid,nstat,lat,long,county endif if(nstate.gt.istate) go to 30 go to 20 30 close (11) if(iscnt.gt.0) go to 50 40 write(*,720) read(*,660) yc if(yc.eq.'y'.or.yc.eq.'Y') go to 10 go to 230 50 write(*,730) cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX c rjq 51 continue cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX iname = istate if(iname.eq.66) iname=52 if(iname.eq.91) iname=53 cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX c rjq if (interactive.or.(istatno.eq.(-1))) then go to 52 endif if(inputfile.eq.'XXX') then go to 53 endif open(10,file=inputfile,status='old',err=55) rewind(10) go to 60 cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 52 read(*,440)istatno 53 open(10,file=dr(iname),status='old',err=55) rewind (10) go to 60 55 write(*,*)' Mother Nature does not want to run state ',dr(iname) write(*,*)' Parameters for this state are not loaded on' write(*,*)' this computer. Check directory for file ',dr(iname) write(*,*)' -- Enter q to quit or c to Continue' read(*,660) yc if(yc.eq.'q'.or.yc.eq.'Q') then go to 230 else go to 10 endif c Gettin Station from State File 60 read(10,450,end=40)stidd,nst,nstat cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX c rjq c if(nst.ne.istate) then if(istate.ne.(-1).and.nst.ne.istate) then read(10,460) go to 60 endif if(nstat.ne.istatno) then read(10,460) go to 60 endif c Read Pricipitation, Temperature, Radiation, etc. c ylt - Lattitude. c yll - Longitude. c rst(i,j) - Array Of Monthly precipitation stats. c obmx - Maximum Temperature. c obmn - Minimum Temperature. c stdmx - Standard Dev. Tmax. c stdmn - Standard Dev. Tmin. c stdsl - Standard Dev. Sol. c wi - Average Maximum .5 Hour Precip. c rh - Dew Point Temp. read(10,470)ylt,yll,itype,tp6 read(10,480)(rst(i,1),i=1,12),(rst(i,2),i=1,12),(rst(i,3),i=1,12) read(10,480)(prw(1,i),i=1,12),(prw(2,i),i=1,12) read(10,480)(obmx(i),i=1,12) read(10,480)(obmn(i),i=1,12) read(10,480)(stdtx(i),i=1,12),(stdtm(i),i=1,12) read(10,480)(obsl(i),i=1,12) read(10,480)(stdsl(i),i=1,12) read(10,480)(wi(i),i=1,12) read(10,480)(rh(i),i=1,12) cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX c John Tatarko 07/16/93 sumtmp = 0. do 1240 i = 1, 12 sumtmp = sumtmp + ( obmx(i) + obmn(i)) c XXXXXXX jt 02/03/98 motemp(i) = (obmx(i) + obmn(i)) / 2. motemp(i) = (motemp(i) - 32.) / 1.8 c XXXXXXX jt end 02/03/98 1240 continue mairtmp = ((sumtmp/24.) - 32.) / 1.8 cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX c Wind Data read in here. c wvl(i,j,k) - array of wind paramters where c i - ith direction (1 - north - 16 nnw) c j - parameters (1 - 4) c 1 - % time from direction i c 2 - mean speed from direction i c 3 - standard deviation of speed from direction i c 4 - skew coeficient of speed from direction i c Calm is treated seperately as direction 0, speed 0 c Only a % time values is need for calm generation. read(10,1250)(((wvl(i,j,k),k=1,12),j=1,4),i=1,16) 1250 format(8x,12f6.2) read(10,1250)(calm(i),i=1,12) read(10,1260)site(1),wgt(1),site(2),wgt(2),site(3),wgt(3) 1260 format(a19,f6.3,2(2x,a19,f6.3)) wmax=0.0 do 70 i=1,12 wi(i)=wi(i) cvtx(i)=stdtx(i)/obmx(i) cvtm(i)=stdtm(i)/obmn(i) if(obsl(i).le.0.0) then cvs(i)=0.0 go to 70 endif cvs(i)=stdsl(i)/obsl(i) 70 continue do 80 i=1,12 wi(i)=wi(i) 80 continue do 6000 i=1,12 6000 dir(i,1)=wvl(1,1,i) do 7010 i=1,12 do 7000 j=2,16 dir(i,j) = dir(i,j-1)+wvl(j,1,i) 7000 continue j=17 dir(i,j)=100.0 7010 continue do 7050 i=1,12 do 7050 j=1,17 dir(i,j)=dir(i,j)*.01 7050 continue close (10) c Close the Parameter File and Write the Parameter if He Wants. cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX c rjq if (.not.interactive) go to 110 cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX write(*,490)stidd,nstat read(*,660)yc if(yc.eq.'n'.or.yc.eq.'N') go to 110 write(*,500)obmx,obmn write(*,520)obsl write(*,540) do 90 i=1,12 write(*,560) prw(1,i),prw(2,i) 90 continue write(*,570) do 100 i=1,12 write(*,600)(rst(i,j),j=1,3) 100 continue write(*,610) write(*,620)(stdtx(i),i=1,12) write(*,620)(stdtm(i),i=1,12) write(*,630)(stdsl(i),i=1,12) write(*,590) write(*,620)(cvtx(i),i=1,12) write(*,620)(cvtm(i),i=1,12) write(*,620)(cvs(i),i=1,12) write(*,620)(wi(i),i=1,12) write(*,265)rh write(*,620)tp6 write(*,2070)site(1),wgt(1),site(2),wgt(2),site(3),wgt(3) 2070 format(1x,'Wind Data Interpolated from', 1 /1x,a19,f6.3,2(2x,a19,f6.3)) write(*,2060)(((wvl(i,j,k),k=1,12),j=1,4),i=1,16) 2060 format(1x,12f6.2) write(*,2060)(calm(i),i=1,12) c c Begin Climate Generation c 110 sml=0.0 r5max=0.0 do 120 i=1,12 if(wi(i).ge.r5max) r5max=wi(i) 120 continue do 125 i=1,12 wi(i)=wi(i) 125 continue call r5mon(tp6) ab=0.02083 ab1=1.0-ab nt=0 clt=57.296 pit=58.13 pi2=6.283185 xx=ylt/clt yls=sin(xx) ylc=cos(xx) vv=randn(k1) l=2 if(vv.gt.prw(1,1)) l=1 rn1=randn(k7) v1=randn(k2) v3=randn(k3) v5=randn(k4) v7=randn(k5) v9=randn(k8) v11=randn(k9) msim=1 nsim=1 cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX c Rob Quinn c Assumption: ibyear, numyear, and outfil were all c specified on the 'command lines'! if (.not.interactive) then iopt=5 iyear=ibyear go to 137 endif cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX c Get Options from User 56 write(*,750) read(*,*)iopt if (iopt.eq.8) go to 230 if (iopt.eq.6) then write(*,*)'Enter Observed Data Input File Name' read(*,640)infile open(9,file=infile,status='old',err=57) rewind(9) read(9,775)ioyr backspace (9) endif go to 58 57 write(*,*)' Error - File not Found ' write(*,*)' Check Directory for Observed Data File or' write(*,*)' Renter File Name' write(*,*)' Enter c to continue or q to quit' read(*,'(a1)')yc if(yc.eq.'q'.or.yc.eq.'Q') go to 230 go to 56 58 numyr=1 if (iopt.eq.1) go to 140 c Single storm input data section if (iopt.eq.4.or.iopt.eq.7) then write(*,*) write(*,*)' Enter Month Day and Year of Storm (mo da yr)' read(*,*)mo,jd,ibyear write(*,*)' Enter Return Period in Years ( e.g. 10 for 10 yr.)' read(*,*)rtyear write(*,*)' Storm Duration in Hours (e.g. 6)' read(*,*)rdur if(iopt.eq.7) then write(*,*)' Enter Design Storm Amount (e.g. 6.30)' read(*,*)damt endif px=1.0 - 1.0/rtyear pxt=px call nrmd(pxt) go to 130 endif if(iopt.eq.6) then ibyear=ioyr numyr=100 go to 130 endif write(*,755) read(*,*) ibyear write(*,760) read(*,*) numyr 130 iyear=ibyear if (iopt.eq.2) go to 140 135 write(*,770)istatno read(*,640)outfil c Ouput Options c 3 - CREAMS - GLEAMS c 4 - Single Storm WEPP c 5 - Continuous Storms WEPP c 6 - Observed Data WEPP c 7 - Design Storm WEPP c cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX c Rob Quinn 137 continue cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX if(iopt.eq.3) then open(8,file=outfil,status='new',err=138) rewind (8) go to 140 else open(7,file=outfil,status='unknown',err=138) rewind (7) go to 140 endif 138 write(*,*)' ' write(*,*)' **** File Already Exists *****' write(*,*)' ' write(*,*)' ',outfil write(*,*)' ' write(*,*)' Enter New File Name ' write(*,*)' ' go to 135 c Get Everything Ready to Start Generation by Options Selected 140 do 150 i=1,12 sumpp(i)=0.0 sumptx(i)=0.0 sumptm(i)=0.0 sumprd(i)=0.0 sumpdr(i)=0.0 150 continue sumpp(13)=0.0 c See What Option was Selected and Set Paths nbt=1 if(iopt .ge. 4)then isim=1 if(iopt.eq.4.or.iopt.eq.7) isim =2 cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX c rjq c write(7,778)isim if(interactive.or.doheader) write(7,778)isim c WEPP MODEL NOT READY FOR THIS YET (3 Lines of Header Removed ADN) c LINES REMOVED FOR OPTIONS 4 AND 7 ONLY (9/11/89 ADN) cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX c rjq c if(iopt.ge.4) then if(iopt.ge.4.and.(interactive.or.doheader)) then cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX c John Tatarko 07/16/93 c write(7,644)stidd,ylt,yll write(7,644)stidd,ylt,yll,mairtmp c John Tatarko 02/02/98 write(7,645) (motemp(i), i = 1,12) cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX write(7,646) iyear,numyr write(7,648) endif endif if(iopt.eq.4.or.iopt.eq.7) then nt=0 if((iyear-iyear/4*4).eq.0) nt=1 ntd1 = jdt(nc,jd,mo,nt) nbt = ntd1 numyr=1 endif c Generate Data by Number of Years and Option - MAIN LOOP do 220 ii=1,numyr 160 ntd=365 if(iopt.le.3.or.iopt.eq.5.or.iopt.eq.6) then cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX c rjq c if((iyear-iyear/4*4).eq.0) ntd=366 if((iyear-iyear/4*4).eq.0.and.(.not. 1 ((iyear-iyear/100*100).eq.0.and.(.not. 1 (iyear-iyear/400*400).eq.0)))) ntd=366 c Leap year if divisible by 4 except when (divisible by 100 except c when divisible by 400). (div4 & !(div100 & !div400)) c I think this is the only time option 5 checks for leap years, c so this is the only place I changed it. You might want to change c other options to work with it fully too. c rjq endif do 170 i=1,12 do 170 jk=1,31 prcip(i,jk)=0.0 tgmx(i,jk)=0.0 tgmn(i,jk)=0.0 radg(i,jk)=0.0 dur(i,jk)=0.0 170 continue c Start of the Daily Generation Loop with nbt and ntd if(iopt.eq.4.or.iopt.eq.7) ntd=ntd1 do 180 ida=nbt,ntd if(iopt.eq.6) then msim=0 nsim=0 read(9,777,end=225)irida,itmxg,itmng if(irida.eq.9999) nsim=1 if(itmxg.eq.9999) msim=1 if(itmng.eq.9999) msim=1 r(ida)=irida*.01 tmxg=itmxg tmng=itmng endif idr=ida call jlt(mo,jd,idr,ntd) call clgen call windg th=th*clt wind(ida)=wv windd(ida)=th prcip(mo,jd)=r(ida) tgmx(mo,jd)=tmxg tgmn(mo,jd)=tmng tmxg=(tmxg-32.0)*(5.0/9.0) tmng=(tmng-32.0)*(5.0/9.0) tdp = (tdp-32.0)*(5.0/9.0) radg(mo,jd)=ra if(r(ida).le.0.0) then r(ida)=0.0 dur(mo,jd)=0.0 else call alph dur(mo,jd)=4.605/(-2.0*alog(1.0-r1)) if(dur(mo,jd).gt.24.0) dur(mo,jd)=24.0 endif c Set duration if a single storm is selected if(iopt.eq.4.or.iopt.eq.7) dur(mo,jd)=rdur if (iopt.ge.4) then if(r(ida).gt.0.) then call alph xr=r(ida)*25.4 tpr=0.4*(1.0-exp(-.3/r1))/dur(mo,jd) if(tpr.gt.1.0) tpr =0.99 r5u=(1.0 -exp(-125./xr))*xr r5p=-2.0*xr*alog(1.0-r1) if(r5p.gt.tymax(itype)) r5p = tymax(itype) xmav=r5p/(xr/dur(mo,jd)) if((tmxg+tmng)/2.0.le.0.0) xmav = 1.01 if(xmav.lt.1.01) xmav =1.01 else xr=r(ida)*25.4 xmav=0.0 tpr=0.0 endif if(iopt.eq.4.or.iopt.eq.7) then xmav=dmxi(itype)*rdur/24. if(xmav.lt.1.01) xmav=1.01 tpr=dtp(itype) endif if(iopt.eq.7) then xr = damt*25.4 endif c Write WEPP Continuous Storm File write(7,650) jd,mo,iyear,xr,dur(mo,jd),tpr, 1 xmav,tmxg,tmng,radg(mo,jd),wv,th,tdp endif 180 continue c End of Daily Loop if(iopt.gt.1) go to 190 c Option 1 Stuff write(*,780) read(*,660) yc if(yc.eq.'y'.or.yc.eq.'Y') then iview=1 call clmout(iview) end if write(*,790) read(*,660)yc if(yc.eq.'y'.or.yc.eq.'Y') go to 160 write(*,800) read(*,660) yc if(yc.eq.'y'.or.yc.eq.'Y') then go to 10 endif go to 230 190 continue if(iopt.eq.3) go to 205 if(iopt.ge.4) go to 215 c Option 2 Stuff if(iopt.eq.2) then jj=ii call clmout(0) write(*,820) iyear write(*,*)' ' write(*,670) write(*,240)jj,(sump(i),i=1,12) write(*,250)jj,(sumtx(i),i=1,12) write(*,260)jj,(sumtm(i),i=1,12) write(*,270)jj,(sumrd(i),i=1,12) write(*,280)jj,(sumdr(i),i=1,12) write(*,300)sump(13),stidd 195 an=numyr do 200 i=1,12 sumpp(i)=sumpp(i)+sump(i)/an sumptx(i)=sumptx(i)+sumtx(i)/an sumptm(i)=sumptm(i)+sumtm(i)/an sumprd(i)=sumprd(i)+sumrd(i)/an sumpdr(i)=sumpdr(i)+sumdr(i)/an 200 continue sumpp(13)=sumpp(13)+sump(13)/an endif go to 215 205 continue c Option 3 Stuff if(iopt.eq.3) then nb=1 ne=10 do 210 j=1, 36 write(8,680)nstat,iyear,(r(i),i=nb,ne),j nb=ne + 1 ne=nb + 9 210 continue nb=361 ne=366 if(((iyear/4*4)-iyear).ne.0) r(366)=0.0 j=37 write(8,690)nstat,iyear,(r(i),i=nb,ne),j endif 215 iyear=iyear+1 220 continue c End - MAIN LOOP if(iopt.eq.2) then write(*,290)numyr write(*,670) write(*,310)(sumpp(i),i=1,12) write(*,320)(sumptx(i),i=1,12) write(*,330)(sumptm(i),i=1,12) write(*,340)(sumprd(i),i=1,12) write(*,350)(sumpdr(i),i=1,12) write(*,305)numyr,sumpp(13),stidd endif c Check for Another Run or End cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX c Rob Quinn 225 if (.not.interactive) then stop endif cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX write(*,810) read(*,660) yc if(yc.eq.'y'.or.yc.eq.'Y') go to 10 if(iopt.ge.4) then write(7,*) ' ' close (7) endif if(iopt.eq.3) then write(8,*)' ' close (8) endif 230 stop 'Program terminated by user - returning to DOS' 240 format(1x,'prcp',i3,12f6.2) 250 format(1x,'tmax',i3,12f6.2) 260 format(1x,'tmin',i3,12f6.2) 265 format(1x,'Average Monthly Dew Point Temperature',/1x,12f6.2) 270 format(1x,'rad ',i3,12f6.1) 280 format(1x,'dur ',i3,12f6.2) 290 format(/22x,'Average Values for ',i2,' Years'/) 300 format(/1x,'Annual Precipitation =',f6.2,a30/) 305 format(/1x,'Average Annual Precipitation for ',i2, 1' Years =',f6.2,a30/) 310 format(1x,'prcp ',12f6.2) 320 format(1x,'tmax ',12f6.2) 330 format(1x,'tmin ',12f6.2) 340 format(1x,'rad ',12f6.1) 350 format(1x,'dur ',12f6.2) 360 format(1x,'01 Alabama 20 Michigan 39 S. Dakota', 1/1x,'02 Arizona 21 Minnesota 40 Tennessee', 2/1x,'03 Arkansas 22 Mississippi 41 Texas', 3/1x,'04 California 23 Missouri 42 Utah', 4/1x,'05 Colorado 24 Montana 43 Vermont', 5/1x,'06 Connecticut 25 Nebraska 44 Virginia', 6/1x,'07 Delaware 26 Nevada 45 Washington', 7/1x,'08 Florida 27 New Hampshire 46 West Virginia', 8/1x,'09 Georgia 28 New Jersey 47 Wisconsin', 9/1x,'10 Idaho 29 New Mexico 48 Wyoming') 370 format(1x,'11 Illinois 30 New York 49 Washington DC', 1/1x,'12 Indiana 31 North Carolina 50 Alaska', 2/1x,'13 Iowa 32 North Dakota 51 Hawaii', 3/1x,'14 Kansas 33 Ohio 66 Puerto Rico', 4/1x,'15 Kentucky 34 Oklahoma 91 Pacific Islands', 5/1x,'16 Lousiana 35 Oregon ', 6/1x,'17 Maine 36 Pennsylvania ', 7/1x,'18 Maryland 37 Rhode Island ', 8/1x,'19 Massachusetts 38 South Carolina ') 380 format(1x,'State Climate Code - Available Stations'/) 390 format(i2) 400 format(a32) 410 format(a41,i2) 420 format(a41,i2,i4,i3,i4,a20) 430 format(1x,a41,2x,i4,3x,i3,2x,i4,1x,a19) 440 format(i4) 450 format(a41,i2,i4) 460 format(////////////////////////////////////////////////// 1 //////////////////////////////) 470 format(6x,f6.2,6x,f6.2,16x,i2/33x,f4.2) 480 format(8x,12f6.2) 490 format(/1x,'Do you want to view data found for station',//, 1 a41,2x,i4,' (y/n)?: ') 500 format(/1x,'observed monthly ave max temperature',//, 1 1x,12f6.2,//, 1 1x,'observed monthly ave min temperature',//, 1 1x,12f6.2/) 510 format(12f6.1) 520 format(/1x,'observed monthly ave solar radiation',//, 1 1x,12f6.1/) 530 format(20i4) 540 format(/1x,'wet-dry state probabilities'/) 550 format(4f10.5) 560 format(1x,2f10.5) 570 format(/1x,'mean,st.dev.,and skew coef. of daily rainfall'/) 580 format(3f10.5) 590 format(/1x,'coefficient of variation for max, min temp,', 1 ' solar radiation, and max .5 hr rain'/) 600 format(1x,3f10.5) 610 format(/1x,'standard deviation for max and min temp,', 1 ' and solar radiation'/) 620 format(1x,12f6.2) 630 format(1x,12f6.1) 640 format(a32) cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX c John Tatarko 07/16/93 c 644 format(' Station: ',a41,6x,' CLIGEN VERSION 2.1',/ c 1 ' Latitude Longitude',/,2f9.2) 644 format(' Station: ',a41,6x,' CLIGEN VERSION 2.1a',/ 1 ' Latitude Longitude Avg Annual Air Temp (c)',/, 2 2f9.2, 14x, f7.2) 645 format (12(1x,f6.2)) cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX c rjq c646 format(' Beginning year: ',i3, c 1 ' Number of years simulated: ',i3) 646 format(' Beginning year: ',i4, 1 ' Number of years simulated: ',i4) cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX c rjq c648 format(' day mon year prcp dur tp ip tmax', 648 format(' day mon year prcp dur tp ip tmax', 1 ' tmin rad w-vel w-dir dew', 2 /,16x,'(mm) (h)',15x,'(C) (C)', 3 ' (ly/day) m/sec deg (C)') cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX c rjq c Make the format except 4 digit years c650 format(1x,3i4,1x,2f6.2,f5.2,1x,f6.2,2f7.2,f7.2,f6.2,2f7.2) 650 format(2x,i2,2x,i2,1x,i4,1x,2f6.2,f5.2,1x,f6.2,2f7.2, 1 f7.2,f6.2,2f7.2) cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 660 format(a1) 670 format(1x,'elem',' yr',' J F M A M J', 1 ' J A S O N D'/) 680 format(i5,i5,10f5.2,8x,i2) 690 format(i5,i5,6f5.2,28x,i2) 700 format(/1x,'Enter state climate code (ex. 01 for Alabama): ') 710 format(/1x,'Stations Available',19x,' Station No. ', 1 'Lat. Long County'/, 11x,8('-'),1x,9('-'),19x,17('-'),2x,4('-'),1x,7('-')/) 720 format(/1x,'No stations available - do you want to' 1' continue (y/n)? ') 730 format(/1x,'Enter the station index: ') 740 format(/1x,'Enter (y/n): ') 750 format(/1x,'Weather Generator Options',/,1x,7('-'),1x, 19('-'),1x,7('-'),//,1x,'1 - Single Year Simulation - Screen',/, 11x,'2 - Multiple Year - Screen Output',/, 11x,'3 - Multiple Year Simulation - CREAMS - GLEAMS Output File',/, 11x,'4 - Selected Single Storm WEPP - Output File',/, 11x,'5 - Multiple Year - WEPP Output File',/, 11x,'6 - Read Observed P and Temp and Generate Missing Data',/, 11x,'7 - Single Design Storm - TR 55 Storm Type WEPP Output File',/, 11x,'8 - Exit Weather Generator Program',//, 11x,'Enter generator option (1-8): ') 755 format(/1x,'Enter beginning simulation year', 1' (last 2 digits, ex. 80): ') 760 format(/1x,'Enter number of years to simulate: ') c 1' (last 2 digits, ex. 88): ') 765 format(/1x,'Number of years should be .le. 20' 1' - reenter') 770 format(/1x,'Enter output file name (ex. Indy.dat): ', 1 'Station No. ',i4) 775 format(13x,i2) 777 format(15x,3i5) 778 format(i4) 780 format(/1x,'Do you want to view generated data (y/n)? ') 790 format(/1x,'Do you want to simulate another year (y/n)? ') 800 format(/1x,'Do you want another Station (y/n)? ') 810 format(/1x,'Do you want to continue (y/n)? ') 820 format(/15x,'Summary of Elements Generated - Year 19',i2) end c c c subroutine alph c c This subroutine computes alpha, a dimensionless parameter that c expresses the fraction of total rainfall that occurs c during 0.5 ho. c c common /bk1/ wvl(16,4,12),dir(12,17),wv,th,pi2,ang,rh(12), 1 tdp,calm(12) common /bk3/ alai(5),salb(5),dm(5),rd(5),uw,ep,es,rain,ws, 1 je,j,ida common /bk4/ nc(13),iyr,nt,mo,px,iopt,dtp(4),dmxi(4) common /bk5/ r(366),vmx(5),v(5),flu(6),bp(5),hc(5),hcr(5),fp(5), 1 cs(5),sk(5),smx(5),ev,sp,rl,q1,o,qd,af,pq,yp,aff, 1 amps,cn,sml,cfp(5),swm(5),sfc(5),dk(5) common /bk7/ rst(12,3),t0(5),prw(2,12),obmx(12),obmn(12), 1 obsl(12),wft(12),subp(5),sno(5),snoev,k1(4), 1 k2(4),cvs(12),k3(4),k4(4),k5(4),k6(4),k7(4), 1 k8(4),k9(4),cvtx(12),cvtm(12),v1,v3,v5,v7,v9,v11,ra, 1 tmxg,tmng,rmx,yls,ylc,pit,st0,amp,avt,ts, 1 nsim,msim,l,tmx(5),tmn(5),tx(5),tb(5),to(5), 1 dmm(12,5),dm1(5),stdtx(12),stdtm(12),stdsl(12) common /bk9/ hyd(250,5),xi(2,5),wi(12),tc(6),al(5),ab,ab1,rn1,r1, 1 pr,lu1,vo,nmb,tp,tr,ddt,tt c c dimension k7(4) c ei=r(ida)-sml ai=ab1/(wi(mo)-ab) if (ei.lt.1.0) go to 10 ajp=1.0-exp(-5.0/ei) go to 20 10 ajp=1.0 20 r1=dstg(rn1,ai,k7) r1=(ei*(ab+r1*(ajp-ab))+sml*ab)/r(ida) return end c c c block data c c contains generator seeds for the weather generator c common /bk1/ wvl(16,4,12),dir(12,17),wv,th,pi2,ang,rh(12), 1 tdp,calm(12) common /bk3/ alai(5),salb(5),dm(5),rd(5),uw,ep,es,rain,ws, 1 je,j,ida common /bk4/ nc(13),iyr,nt,mo,px,iopt,dtp(4),dmxi(4) common /bk7/ rst(12,3),t0(5),prw(2,12),obmx(12),obmn(12), 1 obsl(12),wft(12),subp(5),sno(5),snoev,k1(4), 1 k2(4),cvs(12),k3(4),k4(4),k5(4),k6(4),k7(4), 1 k8(4),k9(4),cvtx(12),cvtm(12),v1,v3,v5,v7,v9,v11,ra, 1 tmxg,tmng,rmx,yls,ylc,pit,st0,amp,avt,ts, 1 nsim,msim,l,tmx(5),tmn(5),tx(5),tb(5),to(5), 1 dmm(12,5),dm1(5),stdtx(12),stdtm(12),stdsl(12) common /bk5/ r(366),vmx(5),v(5),flu(6),bp(5),hc(5),hcr(5),fp(5), 1 cs(5),sk(5),smx(5),ev,sp,rl,q1,o,qd,af,pq,yp,aff, 1 amps,cn,sml,cfp(5),swm(5),sfc(5),dk(5) common /bk9/ hyd(250,5),xi(2,5),wi(12),tc(6),al(5),ab,ab1,rn1,r1, 1 pr,lu1,vo,nmb,tp,tr,ddt,tt common /cl1/ prcip(12,31),tgmx(12,31),tgmn(12,31),radg(12,31), 1 dur(12,31) common /sumr/sump(13),sumtx(12),sumtm(12),sumrd(12), 1 sumdr(12) c data k1/9,98,915,92/ data k2/135,28,203,85/ data k3/43,54,619,33/ data k4/645,9,948,65/ data k5/885,41,696,62/ data k6/51,78,648,0/ data k7/227,57,929,37/ data k8/205,90,215,31/ data k9/320,73,631,49/ data nc/0,31,59,90,120,151,181,212,243,273,304,334,365/ data dtp/.4,.32,.5,.5/ data dmxi/18.24,5.76,32.88,20.16/ end c c c subroutine clgen c c This subroutine simulates daily solar radiation, simulates daily c precipitation and/or maximum and minimum air temperature at the c users option, and calls functions randn and dstn1. c common /bk1/ wvl(16,4,12),dir(12,17),wv,th,pi2,ang,rh(12), 1 tdp,calm(12) common /bk3/ alai(5),salb(5),dm(5),rd(5),uw,ep,es,rain,ws, 1 je,j,ida common /bk4/ nc(13),iyr,nt,mo,px,iopt,dtp(4),dmxi(4) common /bk7/ rst(12,3),t0(5),prw(2,12),obmx(12),obmn(12), 1 obsl(12),wft(12),subp(5),sno(5),snoev,k1(4), 1 k2(4),cvs(12),k3(4),k4(4),k5(4),k6(4),k7(4), 1 k8(4),k9(4),cvtx(12),cvtm(12),v1,v3,v5,v7,v9,v11,ra, 1 tmxg,tmng,rmx,yls,ylc,pit,st0,amp,avt,ts, 1 nsim,msim,l,tmx(5),tmn(5),tx(5),tb(5),to(5), 1 dmm(12,5),dm1(5),stdtx(12),stdtm(12),stdsl(12) common /bk5/ r(366),vmx(5),v(5),flu(6),bp(5),hc(5),hcr(5),fp(5), 1 cs(5),sk(5),smx(5),ev,sp,rl,q1,o,qd,af,pq,yp,aff, 1 amps,cn,sml,cfp(5),swm(5),sfc(5),dk(5) c xi=ida sd=0.4102*sin((xi-80.25)/pit) ch=-yls*tan(sd)/ylc if (ch.ge.1.0) go to 10 if (ch.le.-1.0) go to 20 h=acos(ch) go to 30 10 h=0.0 go to 30 20 h=3.1416 30 dd=1.0+0.0335*sin((xi+88.2)/pit) ys=yls*sin(sd) yc=ylc*cos(sd) rmx=711.0*(h*ys+yc*sin(h)) sum=prw(1,mo)*(1.0-prw(1,mo))+prw(2,mo)*(1.0-prw(2,mo)) if(sum.eq.0.0) sum = 0.10 if (nsim.eq.0) go to 40 vv=randn(k1) if (prw(l,mo).le.0.0) go to 60 if (vv.gt.prw(l,mo)) go to 60 v8=randn(k5) if(rst(mo,3).eq.0.0) rst(mo,3)=0.01 r6=rst(mo,3)/6.0 xlv=(dstn1(v7,v8)-r6)*r6+1.0 c Selected storm frequency here if(iopt.eq.4) then call nrmd(px) xlv=(px-r6)*r6+1.0 endif c * * * * * * * * * * * * * * * * xlv=(xlv**3-1.0)*2.0/rst(mo,3) v7=v8 r(ida)=xlv*rst(mo,2)+rst(mo,1) if (r(ida).lt.0.01) r(ida)=0.01 go to 50 40 if (r(ida).le.0.0) go to 60 50 xx=1.0-prw(l,mo) l=1 go to 70 60 xx=prw(l,mo) r(ida)=0.0 l=2 70 if(msim.eq.0) then xx=xx/sum v12=randn(k9) tdp=xx*dstn1(v11,v12) v11=v12 go to 80 endif xx=xx/sum if(xx.gt.1.0) xx = 1.0 v2=randn(k2) tmxg=xx*dstn1(v1,v2) v1=v2 v4=randn(k3) tmng=xx*dstn1(v3,v4) v12=randn(k9) tdp =xx*dstn1(v11,v12) v3=v4 v11=v12 tmxg=obmx(mo)*(1.0+tmxg*cvtx(mo)) tmng=obmn(mo)*(1.0+tmng*cvtm(mo)) if (tmng.gt.tmxg) tmng=tmxg-0.2*abs(tmxg) c Change made to tdp average cvtm and cvtx ADN 80 tdp =rh(mo)*(1.0+tdp*(cvtm(mo)+cvtx(mo)/2.)) if (tdp.gt.((tmxg+tmng)/2.)) tdp=((tmxg+tmng)/2.)*0.99 c Change to limit - dew point in cold months ADNe if(tdp.lt.-10.) tdp=1.1*tmng v6=randn(k4) ra=xx*dstn1(v5,v6) rx=rmx-obsl(mo) if (obsl(mo).gt.rx) rx=obsl(mo) ra=obsl(mo)+ra*rx/4.0 if(ra.ge.rmx) ra=0.90*rmx if (ra.le.0.0) ra=0.05*rmx v5=v6 return end c c c subroutine clmout(iview) c This subroutine calculates daily and monthly values for options c 1 and 2. It writes ouput to screen for option 1. common /bk1/ wvl(16,4,12),dir(12,17),wv,th,pi2,ang,rh(12), 1 tdp,calm(12) common /bk4/ nc(13),iyr,nt,mo,px,iopt,dtp(4),dmxi(4) common /cl1/ prcip(12,31),tgmx(12,31),tgmn(12,31),radg(12,31), 1 dur(12,31) common /sumr/sump(13),sumtx(12),sumtm(12),sumrd(12), 1 sumdr(12) c do 10 i =1,12 sump(i) = 0.0 sumtx(i)= 0.0 sumtm(i)= 0.0 sumrd(i)= 0.0 sumdr(i)= 0.0 10 continue sump(13)= 0.0 do 20 i =1,12 an = nc(i+1)-nc(i) do 20 j =1,31 sump(i)=sump(i)+prcip(i,j) sump(13)=sump(13)+prcip(i,j) sumtx(i)=sumtx(i)+tgmx(i,j)/an sumtm(i)=sumtm(i)+tgmn(i,j)/an sumrd(i)=sumrd(i)+radg(i,j)/an sumdr(i)=sumdr(i)+dur(i,j) 20 continue if(iopt.eq.2) return if(iview.eq.0) go to 90 30 write(*,850) read(*,*)i goto(40,50,60,70,80,90),i 40 write(*,125) write(*,150) write(*,200)(j,(prcip(i,j),i=1,12),j=1,31) write(*,150) write(*,300)(sump(i),i=1,12) write(*,350)sump(13) go to 30 50 write(*,400) write(*,150) write(*,450)(j,(tgmx(i,j),i=1,12),j=1,31) write(*,150) write(*,500)(sumtx(i),i=1,12) go to 30 60 write(*,550) write(*,150) write(*,600)(j,(tgmn(i,j),i=1,12),j=1,31) write(*,150) write(*,500)(sumtm(i),i=1,12) go to 30 70 write(*,650) write(*,150) write(*,700)(j,(radg(i,j),i=1,12),j=1,31) write(*,150) write(*,500)(sumrd(i),i=1,12) go to 30 80 write(*,750) write(*,150) write(*,800)(j,(dur(i,j),i=1,12),j=1,31) write(*,150) write(*,300)(sumdr(i),i=1,12) go to 30 c 90 continue return c 125 format(20x,'Daily Precipitation Amounts (Inches)'/) 150 format(1x,' J F M A M J J A' 1' S O N D') 200 format(1x,i2,12f6.2) 250 format(/19x,'Single Year Simulation Data Set - Year ',i2/) 300 format(/1x,'tot',f5.2,11f6.2) 350 format(/1x,'Annual Precipitation ',f6.2/) 400 format(/20x,'Daily Maximum Air Temperature (F)'/) 450 format(1x,i2,12f6.0) 500 format(/1x,'ave',f5.0,11f6.0/) 550 format(/20x,'Daily Minimum Air Temperature (F)'/) 600 format(1x,i2,12f6.0) 650 format(/20x,'Daily Solar Radiation (Langleys) '/) 700 format(1x,i2,12f6.0) 750 format(/20x,'Storm Duration (Hours)'/) 800 format(1x,i2,12f6.2) 850 format(/1x,'Climate Data Viewing Options',/,1x,7('-'), 11x,4('-'),1x,7('-'),1x,7('-'),/,1x,'1 - Precipitation',/,1x, 1'2 - Maximum Air Temperature',/,1x,'3 - Minimum Air' 1' Temperature',/,1x,'4 - Solar Radiation',/,1x,'5 - Storm' 1' Duration',/,1x,'6 - End',//,1x,'Enter viewing option') end c c c function dstg (rn1,ai,k7) c c This function provides numbers from a gamma distribution, given c two random numbers. dimension k7(4) c data xn1/10.0/ c 10 dstg=rn1 xx=rn1*ai rn=randn(k7) fu=xx**xn1*exp(xn1*(1.0-xx)) rn1=rn if (fu.lt.rn) go to 10 return end c c c. function dstn1 (rn1,rn2) c c This function computes the distance from the mean of a c normal distribution with mean = 0 and standard deviation c = 1, given two random numbers. c dstn1=sqrt(-2.0*alog(rn1))*cos(6.283185*rn2) c return end c c c function jdt (nc,i,m,nt) c c This subroutine computes the day of the year, given the month and c the day of the month. c dimension nc(13) c if (m.gt.2) go to 10 jdt=nc(m)+i go to 20 10 jdt=nc(m)+nt+i 20 return end c c c subroutine jlt (mo,nday,jday,ntd) c dimension nn(12) data nn/31,28,31,30,31,30,31,31,30,31,30,31/ c if(ntd.eq.366) then nn(2) = 29 else nn(2) = 28 endif kday = jday do 40 j = 1,12 jday = jday -nn(j) if (jday) 30,30,40 30 nday = jday +nn(j) go to 50 40 continue 50 jday = kday mo = j return end c c c subroutine r5mon(tp6) c This subroutine Smothes and corrects .r hour monthly rainfall c intensity data. common /bk9/ hyd(250,5),xi(2,5),wi(12),tc(6),al(5),ab,ab1,rn1,r1, 1 pr,lu1,vo,nmb,tp,tr,ddt,tt common /bk7/ rst(12,3),t0(5),prw(2,12),obmx(12),obmn(12), 1 obsl(12),wft(12),subp(5),sno(5),snoev,k1(4), 1 k2(4),cvs(12),k3(4),k4(4),k5(4),k6(4),k7(4), 1 k8(4),k9(4),cvtx(12),cvtm(12),v1,v3,v5,v7,v9,v11,ra, 1 tmxg,tmng,rmx,yls,ylc,pit,st0,amp,avt,ts, 1 nsim,msim,l,tmx(5),tmn(5),tx(5),tb(5),to(5), 1 dmm(12,5),dm1(5),stdtx(12),stdtm(12),stdsl(12) common /bk4/ nc(13),iyr,nt,mo,px,iopt,dtp(4),dmxi(4) c dimension sm(12),smm(12) c xy2=0.5/tp6 do 10 i=1,12 wi(i)=wi(i)*25.4 10 continue sm(1)=(wi(12)+wi(1)+wi(2))/3.0 do 20 i=2,11 sm(i)=(wi(i-1)+wi(i)+wi(i+1))/3.0 20 continue sm(12)=(wi(11)+wi(12)+wi(1))/3.0 do 30 i=1,12 if(prw(1,i).eq.0.0) prw(1,i)=0.01 if(prw(2,i).eq.0.0) prw(2,i)=0.01 if(rst(i,1).eq.0.0) rst(i,1)=0.01 xm=nc(i+1)-nc(i) smm(i)=xm*prw(1,i)/(1.0-prw(2,i)+prw(1,i)) r25=rst(i,1) f=xy2/smm(i) wi(i)=-sm(i)/alog(f) wi(i)=1.0-exp(-wi(i)/r25) if (wi(i).lt.0.1) wi(i)=0.1 if (wi(i).gt.0.95) wi(i)=0.95 30 continue return end c c c function randn (k) c c This function provides random numbers ranging from 0. to 1.0 c dimension k(4) c k(4)=3*k(4)+k(2) k(3)=3*k(3)+k(1) k(2)=3*k(2) k(1)=3*k(1) i=k(1)/1000 k(1)=k(1)-i*1000 k(2)=k(2)+i i=k(2)/100 k(2)=k(2)-100*i k(3)=k(3)+i i=k(3)/1000 k(3)=k(3)-i*1000 k(4)=k(4)+i i=k(4)/100 k(4)=k(4)-100*i randn=(((float(k(1))*.001+float(k(2)))*.01+float(k(3)))*.001+ 1 float(k(4)))*.01 return end subroutine windg c this subroutine simulates daily average wind velocity using mean, c standard deviation, and skew coefficient of wind speed. c Wind direction is selected from a uniform distribution. common /bk1/ wvl(16,4,12),dir(12,17),wv,th,pi2,ang,rh(12), 1 tdp,calm(12) common /bk3/ alai(5),salb(5),dm(5),rd(5),uw,ep,es,rain,ws, 1 je,j,ida common /bk4/ nc(13),iyr,nt,mo,px,iopt,dtp(4),dmxi(4) common /bk7/ rst(12,3),t0(5),prw(2,12),obmx(12),obmn(12), 1 obsl(12),wft(12),subp(5),sno(5),snoev,k1(4), 1 k2(4),cvs(12),k3(4),k4(4),k5(4),k6(4),k7(4), 1 k8(4),k9(4),cvtx(12),cvtm(12),v1,v3,v5,v7,v9,v11,ra, 1 tmxg,tmng,rmx,yls,ylc,pit,st0,amp,avt,ts, 1 nsim,msim,l,tmx(5),tmn(5),tx(5),tb(5),to(5), 1 dmm(12,5),dm1(5),stdtx(12),stdtm(12),stdsl(12) common /bk5/ r(366),vmx(5),v(5),flu(6),bp(5),hc(5),hcr(5),fp(5), 1 cs(5),sk(5),smx(5),ev,sp,rl,q1,o,qd,af,pq,yp,aff, 1 amps,cn,sml,cfp(5),swm(5),sfc(5),dk(5) c common /bk9/ hyd(250,5),xi(2,5),wi(12),tc(6),al1,ab,ab1,rn1,r1, c 1 pr,lu1,vo,nmb,tp,tr,ddt,tt fx=randn(k6) do 1 j=1,16 j1=j-1 if (dir(mo,j).gt.fx) go to 2 1 continue c calm condition found wv=0.0 th=0.0 return c c Wind direction calculations 2 if (j.eq.1) go to 3 g=(fx-dir(mo,j1))/(dir(mo,j)-dir(mo,j1)) go to 4 3 g=fx/dir(mo,j) 4 xj1=j1 th=pi2*(g+xj1-.5)/16. if (th.lt.0.) th=pi2+th c Wind speed calculations v10=randn(k8) if(wvl(j,4,mo).eq.0.0) wvl(j,4,mo)=0.01 r6=wvl(j,4,mo)/6.0 xlv=(dstn1(v9,v10)-r6)*r6+1.0 xlv=(xlv**3-1.0)*2.0/wvl(j,4,mo) v9=v10 wv=xlv*wvl(j,3,mo)+wvl(j,2,mo) if (wv.lt.0.) wv=.1 c return end subroutine nrmd(r1) c This subroutine returns the standard normal deviate for a c given probability value (e.g. r1=.99 returns r1=2.328) if(r1 .ge. .5) go to 1 sgn = -1. hp = r1 go to 2 1 sgn = 1. hp = 1. - r1 2 t = sqrt(alog(1./(hp*hp))) r1 = sgn*(t-(2.30753+.27061*t)/(1.+.99229*t+.04481*t*t)) return end subroutine header c Subroutine to write header information to screen. write(*,100) 100 format(////,2x,68('*'),/,2x,'*',66x,'*',/,2x,'*',66x,'*', 1 /,2x,'*',14x,'USDA - WATER EROSION PREDICTION', 1 ' PROJECT',13x,'*',/,2x,'*',66x,'*',/,2x,'*',16x, 1 ' WEPP CLIMATE INPUT DATA GENERATOR',16x,'*'/2x,'*', 1 66x,'*',/,2x,'*',20x,'CONTINUOUS SIMULATION', 1 ' AND',21x,'*',/, 1 2x,'*',23x,'SINGLE STORM OPTIONS',23x,'*',/, 1 2x,'*',66x,'*',/, 1 2x,'*',66x,'*',/,2x,'*',27x,'VERSION 2.1',27x,'*',/, 1 2x,'*',66x,'*',/, 1 2x,'*',27x,'OCTOBER 1990',27x,'*',/,2x,'*',66x,'*',/, 1 2x,'*',66x,'*',/,2x,68('*'),//) return end ! This function returns the number of args/options on the commandline ! It is Lahey lf95 PC compiler specific because it uses the getcl() ! subroutine call - LEW function narg() integer narg character*128 cmdline integer i,j i = 0 narg = 1 !include filename - getcl() doesn't call getcl(cmdline) !Lahey compiler specific i = index(cmdline, '-') !get first arg? if (i .ne. 0) then narg = narg + 1 j = i + 1 endif do while (i .ne. 0) i = index(cmdline(j:), ' -') ! get next arg? ! write(6,*) 'i,j,narg is: ', i,j,narg if (i .ne. 0) then narg = narg + 1 j = j+ i + 2 endif end do return end ! This subroutine returns the argument/option specified by argcnt ! from the commandline ! cannot return argcnt = 1 (filename), therefore, we return "" ! This is Lahey lf95 PC compiler specific because it uses the getcl() ! subroutine call - LEW subroutine argopt(argcnt, option) integer argcnt character*128 option integer narg character*128 cmdline integer i,j if (argcnt .eq. 1) then !Can't return filename option = "" goto 100 endif i = 0 narg = 1 !include filename - getcl() doesn't call getcl(cmdline) !Lahey compiler specific i = index(cmdline, '-') !get index to first arg? if (i .ne. 0) then narg = narg + 1 j = i + 1 endif do while (i .ne. 0) i = index(cmdline(j:), ' -') ! get next arg? ! write(6,*) 'i,j,narg is: ', i,j,narg if (i .ne. 0) then if (narg .eq. argcnt) then if (argcnt .eq. 2) then !handle first arg special option = cmdline(j-1:j+i-1) else option = cmdline(j-2:j+i-1) endif goto 100 endif narg = narg + 1 j = j+ i + 2 else !handle last arg/option special if (narg .eq. argcnt) then if (argcnt .eq. 2) then !handle first arg special option = cmdline(j-1:) else option = cmdline(j-2:) endif goto 100 endif endif end do 100 continue ! write(6,*) 'argcnt is: ', argcnt, ' option is: "', option,'"' return end