c ****************************************************************** c * ****** prime version ****** * c * ****** program calparms.f77 ***** * c * Main program to calculate weather generator parameter for the * c * WEPP Project (Water Erosion Predeiction Project * c * Program accesses files created from NWS tapes for each state * c * and calculates the parameters for rainfall and max - min * c * temperature. Daily rainfall and temperature files are on disk * c * in proccessed format. * c * program accesses files, reads data, checks for missing data, * c * summerizes data by month, calculates weather generator * c * parameters then writes completed record to another file for * c * database useage. * c ****************************************************************** c c common /blk1/x(12,31),ss1(12),ssi(12),si(12),sxst(12),nwc(12), 1ndc(12),ntc(12),nrc(12),s1(12),ngn(60),label(20) common /blk2/ sumx(12),ssx(12),sd(12),xmean(12),sk(12),sx3(12) common /blk3/ pmax(12),title(80),nc(12),misda(12),nrec(12) common /blk4/ wc(12),rc(12),dc(12),tc(12),pww(12),pdw(12),pdd(12), 1 pwd(12),xmonth(12),ngrec(12),ntrec(12) common /blk5/ rst(12,3),prw(2,12),k1(4),k2(4),mo,v1,v2,st,lw,nrs common /blkt/ sumt(2,12),sst(2,12),sdt(2,12),xmeant(2,12), 1 nrtc(2,12),misdt(2,12) common /blkx/idx(91,2) common /blksr/ solrad(12),sdsol(12),ax5p(12),dewpt(12),tp5,tp6 dimension nr(366),prc(366) character*32 fil,afil,bfil character*4 el(3),elem character*23 nam character*1 b(32) character*15 county data el/'PRCP','TMAX','TMIN'/ c data k1/9,98,915,92/,k2/135,28,203,85/ nwcode=0 c write(*,*)'Enter climate file name to processed' c read(*,45)fil 45 format(a32) c write(*,*) ' Enter output file name' c read(*,45) afil open (10,file='sindex',status='old') open (11, file='sdata',status='old',recl=69, 1 form='formatted',access='direct') c open (8,file=afil,status='new',access='sequential') c write(*,*) 'Enter Station number' c read(*,46)istat do 5000 k=1,91 read(10,5500)i,(idx(i,j),j=1,2) c write(*,*)k,idx(i,1),idx(i,2) 5000 continue 5500 format(i2,2i8) c46 format(i4) open (6,file='list',status='old',access='sequential',err=2050) 10 read(6,45,end=2060)fil write(*,*)'Reading from file ',fil read(fil,'(32a1)')b b(10)='p' b(11)='a' b(12)='r' write(bfil,'(32a1)')b open (7,file=fil,status='old',access='sequential') open(12,file=bfil,status='new') do 40 i=1,12 sumx(i)=0.0 ssx(i)=0.0 sx3(i)=0.0 pmax(i)=0.0 xmean(i)=0.0 xmonth(i)=0.0 nwc(i)=0 ndc(i)=0 nrc(i)=0 nrec(i)=0 ntrec(i)=0 ntc(i)=0 misda(i)=0 do 40 j=1,2 sumt(j,i)=0.0 sst(j,i)=0.0 sdt(j,i)=0.0 xmeant(j,i)=0.0 misdt(j,i)=0 nrtc(j,i)=0 40 continue years=0.0 tp5=0.0 tp6=0.0 st=0.0 nt=0 l=0 1 read(7,2,end=50)nst,nstat,idist,ny,elem,(nr(i),i=1,366) ist=nst 2 format(i2,i4,i2,i4,a4,366i4) c if(nstat.ne.istat.and.l.eq.0) go to 1 if(elem.ne.el(1)) go to 210 c if(nstat.ne.istat) go to 50 l=1 years=years+1.0 lda = 365 if(((ny/4*4)-ny).eq.0) lda = 366 do 200 i = 1,lda call julin(mo,ida,ny,i) if(nr(i).eq.9999) misda(mo)=misda(mo)+1 ntc(mo)=ntc(mo)+1 if(nr(i).gt.0.and.nr(i).lt.9999) nrc(mo)=nrc(mo)+1 if(nr(i).le.0.and.nt.ne.0) ndc(mo)=ndc(mo)+1 if(nr(i).gt.0.and.nr(i).ne.9999.and.nt.gt.0) 1 nwc(mo)=nwc(mo)+1 nt=nr(i) if(nr(i).eq.9999) nt=0 if(nr(i).eq.0) go to 200 if(nr(i).ge.9998) go to 200 st=nr(i) st=st*.01 x(mo,ida)=st sumx(mo)=sumx(mo)+st ssx(mo)=ssx(mo)+(x(mo,ida)*x(mo,ida)) sx3(mo)=sx3(mo)+(x(mo,ida)*x(mo,ida)*x(mo,ida)) if(x(mo,ida).gt.pmax(mo)) pmax(mo)=x(mo,ida) 200 continue do 201 i=1,12 201 nrec(i)=nrec(i)+1 go to 1 210 ld=1 lda=365 if(((ny/4*4)-ny).eq.0) lda=366 j=1 if(elem.eq.el(3)) j=2 do 220 i=1,lda call julin(mo,ida,ny,i) if(nr(i).eq.9999) misdt(j,mo)=misdt(j,mo) + 1 if(nr(i).lt.9999) nrtc(j,mo)=nrtc(j,mo) + 1 if(nr(i).eq.9999) go to 220 st=nr(i) sumt(j,mo)=sumt(j,mo)+st sst(j,mo)=sst(j,mo)+(st*st) 220 continue go to 1 50 do 250 i=1,12 wc(i)=nwc(i) rc(i)=nrc(i) dc(i)=ndc(i) tc(i)=ntc(i) pww(i)=wc(i)/rc(i) pdw(i)=1.-pww(i) pwd(i)=(rc(i)-wc(i))/(tc(i)-rc(i)) pdd(i)=1.-pdw(i) 250 continue do 300 i=1,12 if(nrc(i).lt.3) nrc(i)=3 an= nrc(i) am=nrec(i) write(*,*)an,am xmonth(i)=sumx(i)/am xmean(i)=sumx(i)/an sd(i)=ssx(i)-(sumx(i)*sumx(i)/an) sd(i)=sqrt(sd(i)/(an-1.0)) sum3=sumx(i)*sumx(i)*sumx(i) sk(i)=an*an*sx3(i)-3.0*an*sumx(i)*ssx(i)+2.0*sum3 sk(i)=sk(i)/(an*(an-1.0)*(an-2.0)) 300 sk(i)=sk(i)/(sd(i)*sd(i)*sd(i)) do 310 j=1,2 do 310 i=1,12 an=nrtc(j,i) xmeant(j,i)=sumt(j,i)/an sdt(j,i)=sst(j,i)-(sumt(j,i)*sumt(j,i)/an) sdt(j,i)=sqrt(sdt(j,i)/(an-1.0)) 310 continue write(*,104)istat,nst 104 format(22x,'station number ',i4,' state 'i2) write(*,100) write(*,101)(i,nrec(i),misda(i),nrc(i),nwc(i),ndc(i),ntc(i), 1 xmonth(i),pmax(i),i=1,12) 100 format(1x,'month num miss. wet w-w dry total mean' 1 /1x,' rec days days days days days rain') 101 format(1x,i4,6i7,2f8.2) write(*,102) 102 format(1x,'month p(w/w) p(w/d) mean std. skew', 1 /1x,' daily dev. coef.') write(*,103)(i,pww(i),pwd(i),xmean(i),sd(i),sk(i),i=1,12) 103 format(1x,i4,f10.3,f8.3,f7.3,f8.3,f8.3) do 130 i=1,12 xmonth(i)=0.0 sumx(i)=0.0 ngrec(i)=0 prw(1,i)=pww(i) prw(2,i)=pwd(i) rst(i,1)=xmean(i) rst(i,2)=sd(i) 130 rst(i,3)=sk(i) write(*,1030) (xmeant(1,i),i=1,12) write(*,1030) (xmeant(2,i),i=1,12) write(*,1030) (sdt(1,i),i=1,12) write(*,1030) (sdt(2,i),i=1,12) write(*,1031) (misdt(1,i),i=1,12) write(*,1031) (misdt(2,i),i=1,12) write(*,1031) (nrtc(1,i),i=1,12) write(*,1031) (nrtc(2,i),i=1,12) 1030 format(1x,12f6.2) 1031 format(1x,12i6) c write parameters to a file call getsta(nst,nstat,nam,mlad,mlam,mlod,mlom,melev, 1 ity,lerr) write(*,*) 'got to here' if(lerr.eq.1) then write(*,*)'station not in station index list' write(*,*)'****** Check Station # ******' nam='station or file err' go to 2060 endif xlat=mlam xlat=xlat/60.0 xlat=mlad+xlat xlon=mlom xlon=xlon/60.0 xlon=mlod+xlon xlon1=-xlon elev=melev*10 c*********************************************************************** c call intsr(xlat,xlon) c c********************************************************************** write(12,4000)nam,nst,nstat,nwcode 4000 format(1x,a23,17x,i2,i4,i2) write(12,4010)xlat,xlon1,years,ity write(12,4020)elev,tp5,tp6 write(12,4030)(xmean(i),i=1,12) write(12,4040)(sd(i),i=1,12) write(12,4050)(sk(i),i=1,12) write(12,4060)(pww(i),i=1,12) write(12,4070)(pwd(i),i=1,12) write(12,4080)(xmeant(1,i),i=1,12) write(12,4085)(xmeant(2,i),i=1,12) write(12,4090)(sdt(1,i),i=1,12) write(12,4095)(sdt(2,i),i=1,12) c********************************************************************* c write(12,4100)(solrad(i),i=1,12) 4100 format(' SOL.RAD',12F6.0) write(12,4110)(sdsol(i),i=1,12) 4110 format(' SD SOL ',12F6.1) write(12,4120)(ax5p(i),i=1,12) 4120 format(' MX .5 P',12F6.2) c********************************************************************* c call intdp(xlat,xlon) c c********************************************************************* write(12,4130)(dewpt(i),i=1,12) 4130 format(' DEW PT ',12F6.2) c c********************************************************************** 4010 format(' LATT=',f7.2,' LONG=',f7.2,' YEARS=',f4.0,' TYPE=',i2) 4020 format(' ELEVATION =',f6.0,' TP5 =',f5.2,' TP6=',f5.2) 4030 format(' MEAN P ',12f6.2) 4040 format(' S DEV P',12f6.2) 4050 format(' SQEW P',12f6.2) 4060 format(' P(W/W) ',12f6.2) 4070 format(' P(W/D) ',12f6.2) 4080 format(' TMAX AV',12f6.2) 4085 format(' TMIN AV',12f6.2) 4090 format(' SD TMAX',12f6.2) 4095 format(' SD TMIN',12f6.2) close (7) close (12) go to 10 2050 write(*,*) 'error opening list file' write(*,*) 'Make sure file list exits in this directory' write(*,*) 'Create with ls *.dat > list ' 2060 close (6) close (7) close (10) c close (11) close (12) stop end subroutine getsta(ist,nstat,nam,mlad,mlam,mlod,mlom,melev, 1 ity,lerr) c fortran program to read tape list inventory files on the prime c and build an interactive file database. c adn 4/29/87 durant, ok common /blkx/idx(91,2) character*23 nam lerr=1 c open (11, file='sdata',status='old',recl=69, c 1 form='formatted',access='direct') c rewind(11) do 560 i=idx(ist,1),idx(ist,2) read(11,570,rec=i)mst,mstat,ndis,nam,mby,mbm,mey,mem,mlad,mlam, 1 mlod,mlom,melev,ity,county 570 format(i2,i4,a1,a23,4i2,2i2,i3,i2,i4,1x,i2,a15) write(*,*)nam,i if(mstat.eq.nstat) then lerr = 0 c close (11) return endif 560 continue 580 format(1x,i2,i8,2x,a23,2x,i2,2x,i2,2x,i3,2x,i2,i8) c close (11) write(*,*)idx(ist,1),idx(ist,2) return end subroutine julin (mo,nday,nyer,jday) dimension ii(12) kday = jday ii(1) = 31 ii(2) = 28 ii(3) = 31 ii(4) = 30 ii(5) = 31 ii(6) = 30 ii(7) = 31 ii(8) = 31 ii(9) = 30 ii(10) = 31 ii(11) = 30 ii(12) = 31 nyt = nyer/ 4 nyt = nyt*4 if ( nyt - nyer ) 20,10,20 10 ii(2) = 29 20 do 30 j = 1,12 jday = jday -ii(j) if (jday) 40,40,30 40 nday = jday +ii(j) go to 35 30 continue 35 jday = kday mo = j return end C This program was written Nov.,1988 by Gene Gander to interpolate values C at specified lat-lon points using closest enclosing triangles of three C points where there exist known values for the parameter of interest. subroutine intsr (tlat,tlon) common /blksr/ solrad(12),sdsol(12),ax5p(12),dewpt(12),tp5,tp6 character*32 sta(15),statin character*5 aa dimension xlat(15),xlon(15),x(15),y(15),far(15),indx(15),ic(3) 1,wgt(3),srx(12),sdx(12),apx(12),dpx(12) 101 format (a32,7x,2f6.2,////) tp5=0.0 tp6=0.0 do 61 i=1,12 solrad(i)=0.0 sdsol(i)=0.0 ax5p(i)=0.0 dewpt(i)=0.0 srx(i)=0.0 sdx(i)=0.0 apx(i)=0.0 61 dpx(i)=0.0 1 do 2 i=1,10 2 far(i)=99999999.0 open(13,file='statparm',status='old') call lltoalbr (tlat,tlon,tx,ty) j=0 3 j=j+1 read(13,101,end=7) statin,alat,alon call distance (tlat,tlon,alat,alon,dist) do 4 i=10,1,-1 if(dist.lt.far(i)) go to 4 k=i+1 go to 5 4 continue k=1 5 do 6 i=10,k,-1 far(i)=far(i-1) sta(i)=sta(i-1) xlat(i)=xlat(i-1) xlon(i)=xlon(i-1) 6 indx(i)=indx(i-1) far(k)=dist xlat(k)=alat xlon(k)=alon sta(k)=statin indx(k)=j go to 3 7 do 8 i=1,10 8 call lltoalbr (xlat(i),xlon(i),x(i),y(i)) ic(1)=0 ic(2)=0 ic(3)=0 do 21 i=3,10 do 21 k=2,i-1 do 21 l=1,k-1 if(x(l)-tx) 10,9,10 9 slpx1=9999. go to 11 10 slpx1=(y(l)-ty)/(x(l)-tx) 11 cptxt1=y(l)-slpx1*x(l) if(x(k)-x(i)) 13,12,13 12 slpki=9999. go to 14 13 slpki=(y(k)-y(i))/(x(k)-x(i)) 14 cptki=y(k)-slpki*x(k) if(slpx1-slpki) 36,35,36 35 sx1=9999.0 go to 37 36 sx1=(cptki-cptxt1)/(slpx1-slpki) 37 sy1=slpx1*sx1+cptxt1 if((sx1.ge.x(k).and.sx1.le.x(i)).or.(sx1.le.x(k).and.sx1.ge.x(i))) 1 then if((sy1.ge.y(k).and.sy1.le.y(i)).or.(sy1.le.y(k).and.sy1.ge.y(i))) 1 then if(x(k)-tx) 16,15,16 15 slpx2=9999. go to 17 16 slpx2=(y(k)-ty)/(x(k)-tx) 17 cptxt2=y(k)-slpx2*x(k) if(x(i)-x(l)) 19,18,19 18 slpil=9999. go to 20 19 slpil=(y(i)-y(l))/(x(i)-x(l)) 20 cptil=y(i)-slpil*x(i) sx2=(cptil-cptxt2)/(slpx2-slpil) sy2=slpx2*sx2+cptxt2 if((sx2.ge.x(i).and.sx2.le.x(l)).or.(sx2.le.x(i).and.sx2.ge.x(l))) 1 then if((sy2.ge.y(i).and.sy2.le.y(l)).or.(sy2.le.y(i).and.sy2.ge.y(l))) 1 then ic(1)=l ic(2)=k ic(3)=i go to 22 endif endif endif endif 21 continue 22 if(ic(1).ne.0) go to 23 dis1=99999999. do 26 k=2,10 do 26 l=1,k-1 if(x(l).eq.x(k)) then slp=9999.0 else slp=(y(l)-y(k))/(x(l)-x(k)) endif 25 antrcp=y(l)-slp*x(l) if(slp) 27,28,27 28 slpnor=9999. go to 30 27 slpnor=-1./slp 30 antnor=ty-slpnor*tx ssx=(antnor-antrcp)/(slp-slpnor) ssy=slpnor*ssx+antnor d=far(k)+far(l)+(((ssx-tx)*(ssx-tx)+(ssy-ty)*(ssy-ty))**0.5)/500. if(d.ge.dis1) go to 26 dis1=d ic(1)=l ic(2)=k 26 continue l=ic(1) k=ic(2) w=far(l)+far(k) wgt(1)=(w-far(l))/w wgt(2)=(w-far(k))/w wgt(3)=0.0 go to 32 23 l=ic(1) k=ic(2) i=ic(3) wgt(1)=1.0/(1.0+far(l)/far(k)+far(l)/far(i)) wgt(2)=wgt(1)*far(l)/far(k) wgt(3)=wgt(1)*far(l)/far(i) 32 write(*,102) indx(l),sta(l),xlat(l),xlon(l),far(l),wgt(1) write(*,102) indx(k),sta(k),xlat(k),xlon(k),far(k),wgt(2) write(*,102) indx(i),sta(i),xlat(i),xlon(i),far(i),wgt(3) 102 format(i4,a22,4f9.3) close(13) open (13,file='statparm',status='old') j=0 51 j=j+1 if (j.ne.indx(l).and.j.ne.indx(k).and.j.ne.indx(i))then read(13,103,end=29) aa go to 51 103 format(a5,////) else read(13,104,end=29) tp5x,tp6x 104 format(59x,2f6.0) read(13,105) (srx(m),m=1,12) read(13,105) (sdx(m),m=1,12) read(13,105) (apx(m),m=1,12) read(13,105) (dpx(m),m=1,12) 105 format (8x,12f6.0) endif if(j.eq.indx(l)) awgt=wgt(1) if(j.eq.indx(k)) awgt=wgt(2) if(j.eq.indx(i)) awgt=wgt(3) do 52 n=1,12 solrad(n)= solrad(n)+awgt*srx(n) sdsol(n)= sdsol(n)+awgt*sdx(n) ax5p(n)= ax5p(n)+awgt*apx(n) 52 dewpt(n)= dewpt(n)+awgt*dpx(n) tp5=tp5+awgt*tp5x tp6=tp6+awgt*tp6x go to 51 29 close(13) return end C This program is just like the one above (subroutine intsr) except c this one was added to handle dewpoint from a different file. subroutine intdp (tlat,tlon) common /blksr/ solrad(12),sdsol(12),ax5p(12),dewpt(12),tp5,tp6 character*32 sta(15),statin character*5 aa dimension xlat(15),xlon(15),x(15),y(15),far(15),indx(15),ic(3) 1,wgt(3),srx(12),sdx(12),apx(12),dpx(12) 101 format (a29,38x,2f3.0,f4.0,f3.0) do 61 i=1,12 dewpt(i)=0.0 61 dpx(i)=0.0 1 do 2 i=1,10 2 far(i)=99999999.0 open(13,file='dewpoint',status='old') call lltoalbr (tlat,tlon,tx,ty) j=0 3 j=j+1 read(13,101,end=7) statin,alatd,alatm,alond,alonm alat=alatd+alatm/60.0 alon=alond+alonm/60.0 call distance (tlat,tlon,alat,alon,dist) do 4 i=10,1,-1 if(dist.lt.far(i)) go to 4 k=i+1 go to 5 4 continue k=1 5 do 6 i=10,k,-1 far(i)=far(i-1) sta(i)=sta(i-1) xlat(i)=xlat(i-1) xlon(i)=xlon(i-1) 6 indx(i)=indx(i-1) far(k)=dist xlat(k)=alat xlon(k)=alon sta(k)=statin indx(k)=j go to 3 7 do 8 i=1,10 8 call lltoalbr (xlat(i),xlon(i),x(i),y(i)) ic(1)=0 ic(2)=0 ic(3)=0 do 21 i=3,10 do 21 k=2,i-1 do 21 l=1,k-1 if(x(l)-tx) 10,9,10 9 slpx1=9999. go to 11 10 slpx1=(y(l)-ty)/(x(l)-tx) 11 cptxt1=y(l)-slpx1*x(l) if(x(k)-x(i)) 13,12,13 12 slpki=9999. go to 14 13 slpki=(y(k)-y(i))/(x(k)-x(i)) 14 cptki=y(k)-slpki*x(k) if(slpx1-slpki) 36,35,36 35 sx1=9999.0 go to 37 36 sx1=(cptki-cptxt1)/(slpx1-slpki) 37 sy1=slpx1*sx1+cptxt1 if((sx1.ge.x(k).and.sx1.le.x(i)).or.(sx1.le.x(k).and.sx1.ge.x(i))) 1 then if((sy1.ge.y(k).and.sy1.le.y(i)).or.(sy1.le.y(k).and.sy1.ge.y(i))) 1 then if(x(k)-tx) 16,15,16 15 slpx2=9999. go to 17 16 slpx2=(y(k)-ty)/(x(k)-tx) 17 cptxt2=y(k)-slpx2*x(k) if(x(i)-x(l)) 19,18,19 18 slpil=9999. go to 20 19 slpil=(y(i)-y(l))/(x(i)-x(l)) 20 cptil=y(i)-slpil*x(i) sx2=(cptil-cptxt2)/(slpx2-slpil) sy2=slpx2*sx2+cptxt2 if((sx2.ge.x(i).and.sx2.le.x(l)).or.(sx2.le.x(i).and.sx2.ge.x(l))) 1 then if((sy2.ge.y(i).and.sy2.le.y(l)).or.(sy2.le.y(i).and.sy2.ge.y(l))) 1 then ic(1)=l ic(2)=k ic(3)=i go to 22 endif endif endif endif 21 continue 22 if(ic(1).ne.0) go to 23 dis1=99999999. do 26 k=2,10 do 26 l=1,k-1 if(x(l).eq.x(k)) then 24 slp=9999. else 25 slp=(y(l)-y(k))/(x(l)-x(k)) endif 33 antrcp=y(l)-slp*x(l) if(slp) 27,28,27 28 slpnor=9999. go to 30 27 slpnor=-1./slp 30 antnor=ty-slpnor*tx ssx=(antnor-antrcp)/(slp-slpnor) ssy=slpnor*ssx+antnor d=far(k)+far(l)+(((ssx-tx)*(ssx-tx)+(ssy-ty)*(ssy-ty))**0.5)/500. if(d.ge.dis1) go to 26 dis1=d ic(1)=l ic(2)=k 26 continue l=ic(1) k=ic(2) w=far(l)+far(k) wgt(1)=(w-far(l))/w wgt(2)=(w-far(k))/w wgt(3)=0.0 go to 32 23 l=ic(1) k=ic(2) i=ic(3) wgt(1)=1.0/(1.0+far(l)/far(k)+far(l)/far(i)) wgt(2)=wgt(1)*far(l)/far(k) wgt(3)=wgt(1)*far(l)/far(i) 32 write(*,102) indx(l),sta(l),xlat(l),xlon(l),far(l),wgt(1) write(*,102) indx(k),sta(k),xlat(k),xlon(k),far(k),wgt(2) write(*,102) indx(i),sta(i),xlat(i),xlon(i),far(i),wgt(3) 102 format(i4,a22,4f9.3) close(13) open (13,file='dewpoint',status='old') j=0 51 j=j+1 if (j.ne.indx(l).and.j.ne.indx(k).and.j.ne.indx(i))then read(13,103,end=29) aa go to 51 103 format(a5) else read(13,105,end=29) (dpx(m),m=1,12) 105 format (31x,12f3.0) endif if(j.eq.indx(l)) awgt=wgt(1) if(j.eq.indx(k)) awgt=wgt(2) if(j.eq.indx(i)) awgt=wgt(3) do 52 n=1,12 52 dewpt(n)= dewpt(n)+awgt*dpx(n) go to 51 29 close(13) return end C This program was written Nov., 88 by Gene Gander. It returns the C distance (dist) between two points of latitude-longitude (degrees) in the C same units as the earth's radius. Presently equal to 6370.997 kilometers. subroutine distance (alat,alon,blat,blon,dist) pi=3.141592654 radius=6370.997 a=alat*pi/180. b=alon*pi/180. d=blon*pi/180. c=blat*pi/180. if(abs(alat).le.abs(blat+.01).and.abs(alat).ge.abs(blat-.01).and. 1abs(alon).le.abs(blon+.01).and.abs(alon).ge.abs(blon-.01)) then dist=0.0 else variab=cos(a)*cos(c)*cos(d-b)+sin(a)*sin(c) if(variab.lt.-1.0) then write(*,*) variab variab=-1.0 endif if(variab.gt.1.0) then write(*,*) variab variab=1.0 endif dist=radius*acos(variab) endif return end C THIS SUBROUTINE CONVERTS LAT-LON POINTS TO X-Y COORDINATES. C RETURNED POINTS ARE FROM AN ALBERS PROJECTION WHERE STANDARD PARALLELS C OF LATITUDE OF 29.5 AND 45.5 ARE USED. THE SUBROUTINE C ASSUMES AN EARTH RADIUS OF 6370997 METERS AND A VERTICAL MERIDIAN C OF LONGITUDE OF 97 DEGREES. IF THESE ASSUMPTIONS ARE NOT VALID THEN C MAKE APPROPRIATE CHANGES. SUBROUTINE LLTOALBR(alat,alon,x,y) REAL LAMDA0,lambda PI=3.141592654 R=6370997. PHI0=90.*PI/180. PHI1=29.5*PI/180. PHI2=45.5*PI/180. LAMDA0=-97.*PI/180. en=(sin(phi1)+sin(phi2))/2. c=cos(phi1)**2+2.*en*sin(phi1) rho0=r*(c-2.*en*sin(phi0))**.5/en lambda=-alon*pi/180. phi=alat*pi/180. theta=en*(lambda-lamda0) rho=r*(c-2.*en*sin(phi))**.5/en x= rho*sin(theta) y=rho0-rho*cos(theta) return END