c ****************************************************************** c * ****** program dat2par ***** * c * (calparms.f77) * c * Main program to calculate weather generator parameter file * c * for the WEPP Project (Water Erosion Prediction 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 * summarizes data by month, calculates weather generator * c * parameters then writes completed record to another file for * c * database useage. * c ****************************************************************** c c modified calparms.for -> dat2par DEH USFS 06/1999 c display non-random third-closest station if only 2 are used. c remove computed gotos from subroutines. c Opens and reads unit 16 file 'list' -- list of .dat files to process c Opens and reads unit 7 file fil -- .dat file c Opens and writes unit 12 file bfil -- .par file c Calls julin c intsr -- triangulate solar radiation plus c intdp -- triangulate monthly dewpoint values from dewpoint file c inttp -- triangulate monthly time-to-peak values c intwind -- triangulate monthly wind data from wind data files c Note: symbol ISTAT referenced but not set c unused symbol COUNTY c unused symbol AFIL c unused symbol PRC !MS$DEBUG common /blk1/x(12,31),ss1(12),ssi(12),si(12),sxst(12),nwc(12), 1 ndc(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),tp(12),wnd(12,4,17),wgt(3) ! ,prc(366) character*32 fil,bfil ! ,afil character*12 stnum character*4 el(3),elem character*23 nam character*1 b(32) c character*15 county character*19 site(3) character*7 point(17),stat(3) character*80 line data el/'PRCP','TMAX','TMIN'/ c data k1/9,98,915,92/,k2/135,28,203,85/ data point/'% N ','% NNE ','% NE ','% ENE ','% E ' 1,'% ESE ','% SE ','% SSE ','% S ','% SSW ','% SW ' 2,'% WSW ','% W ','% WNW ','% NW ','% NNW ','CALM '/ data stat /'MEAN ','STD DEV','SKEW '/ nwcode=0 c write(*,*)'Enter climate file name to be processed' c read(*,45)fil c write(*,*) ' Enter output file name' c read(*,45) afil c open (10,file='sindex',status='old') c open (11, file='sdata',status='old',recl=69, c open (11, file='sdata',status='old',recl=73, c 1 form='formatted',access='direct') c open (8,file=afil,status='new',access='sequential') c write(*,*) 'Enter Station number' c read(*,46)istat c do 5000 k=1,91 c read(10,5500)i,(idx(i,j),j=1,2) c write(*,*)k,idx(i,1),idx(i,2) c 5000 continue c 5500 format(i2,2i8) c46 format(i4) open(16,file='list',status='old',access='sequential',err=2050) 10 read(16,45,end=2060)fil 45 format(a32) print *,'========== Read .dat and create .par ===========' write(*,*)'Reading from file ',fil read(fil,'(32a1)')b b(10)='p' b(11)='a' b(12)='r' write(bfil,'(32a1)')b open (11,file='stations.dat', status='old', access='sequential') 450 read (11,*) stnum ! skip past header records c first line after # is skipped if (stnum(1:1) .eq. '#') goto 450 451 read(11,'(a12,i4,i3,i4,i3,i5,i3,1x,a23)',end=452) * stnum, mlad, mlam, mlod, mlom, melev, itype, nam if (stnum(3:10) .ne. fil(3:10)) goto 451 print *,' Found ',nam print *,' Latitude: ', mlad,':',mlam print *,' Longitude: ', mlod,':',mlom print *,' Elevation: ', melev print *,' itype: ', itype goto 453 452 print *, 'Station not found in stations.dat' close(11) goto 10 453 close(11) c c read(16,'(a)') nam c print *,nam c read(16,*) mlad,mlam,mlod,mlom,melev c open (7,file=fil,status='old',access='sequential',err=454) goto 455 454 print *,' Data file not found for station' goto 10 c open(12,file=bfil,status='new') 455 open(12,file=bfil,status='unknown') 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) 2 format(i2,i4,i2,i4,a4,366i4) ist=nst 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 !Adds day for leap year 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)nwc(mo)=nwc(mo)+1 nt=nr(i) if (nr(i).eq.9999) nt=0 if (nr(i).ne.0 .and. nr(i).lt.9998) then 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) end if 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).ne.9999) then st=nr(i) sumt(j,mo)=sumt(j,mo)+st sst(j,mo)=sst(j,mo)+(st*st) end if 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 c print *,' an,am' do 300 i=1,12 if (nrc(i).lt.3) nrc(i)=3 an = nrc(i) am = nrec(i) c 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 c note, istat not set anywhere (was read from getsta) ********* istat=0 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(' 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 c call getsta(nst,nstat,nam,mlad,mlam,mlod,mlom,melev, c 1 ity,lerr) ity = itype c write(*,*) 'got to here' print *, nst,' ',nstat,' ',nam c if(lerr.eq.1) then c write(*,*)'station not in station index list' c write(*,*)'****** Check Station # ******' c nam='station or file err' c go to 2060 c 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 SOLAR RADIATION DATA and MX .5 P c*********************************************************************** c call intsr up here 'cuz it calculates tp5 & tp6 call intsr(xlat,xlon) write(12,4000)nam,nst,nstat,nwcode 4000 format(1x,a23,17x,i2,i4,i2) write(12,4010)xlat,xlon1,years,ity 4010 format(' LATT=',f7.2,' LONG=',f7.2,' YEARS=',f4.0,' TYPE=',i2) write(12,4020)elev,tp5,tp6 4020 format(' ELEVATION =',f6.0,' TP5 =',f5.2,' TP6=',f5.2) write(12,4030)(xmean(i),i=1,12) 4030 format(' MEAN P ',12f6.2) write(12,4040)(sd(i),i=1,12) 4040 format(' S DEV P',12f6.2) write(12,4050)(sk(i),i=1,12) 4050 format(' SKEW P',12f6.2) write(12,4060)(pww(i),i=1,12) 4060 format(' P(W/W) ',12f6.2) write(12,4070)(pwd(i),i=1,12) 4070 format(' P(W/D) ',12f6.2) write(12,4080)(xmeant(1,i),i=1,12) 4080 format(' TMAX AV',12f6.2) write(12,4085)(xmeant(2,i),i=1,12) 4085 format(' TMIN AV',12f6.2) write(12,4090)(sdt(1,i),i=1,12) 4090 format(' SD TMAX',12f6.2) write(12,4095)(sdt(2,i),i=1,12) 4095 format(' SD TMIN',12f6.2) 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 DEW POINT c********************************************************************* call intdp(xlat,xlon) write(12,4130)(dewpt(i),i=1,12) 4130 format(' DEW PT ',12F6.2) c c********************************************************************* c TIME TO PEAK c********************************************************************* call inttp(xlat,xlon,tp) write(12,4131)(tp(i),i=1,12) 4131 format('Time Pk ',12f6.3) c c********************************************************************** c WIND DATA c********************************************************************** call intwind (xlat,xlon,wnd,wgt,site) do 4134 i=1,16 write(12,4102) point(i),(wnd(j,1,i),j=1,12) do 4134 k=2,4 4134 write(12,4102) stat(k-1),(wnd(j,k,i),j=1,12) write(12,4102) point(17),(wnd(j,1,17),j=1,12) write(12,'(a)') ' ' write(12,'(a)') 'INTERPOLATED DATA (station & weighting factor)' write(12,'(a)') ' ' write(12,'(a)') '---Wind Stations---' write(12,4103) site(1),wgt(1),site(2),wgt(2),site(3),wgt(3) 4102 format(a7,1x,12f6.2) 4103 format(a19,f6.3,2(2x,a19,f6.3)) close (7) c Add station information at end of file. DLS open(22,file="triang.tmp",status="old") write(12,'(a)') '---Solar Radiation and Max .5 P Stations---' read(22,79) line ! sol rad write(12,79) line write(12,'(a)') '---Dewpoint Stations---' read(22,79) line ! dewpt write(12,79) line write(12,'(a)') '---Time Peak Stations---' read(22,79) line ! timepeak write(12,79) line 79 format(a80) close(22) close (12) go to 10 2050 write(*,*) 'error opening list file' write(*,*) 'Make sure file list exists in this directory' write(*,*) 'Create with dir *.dat > list ' 2060 continue c print *,' That''s all' close (6) close (7) close (10) c close (11) close (12) stop end c************************* c GETSTA c************************* subroutine getsta(ist,nstat,nam,mlad,mlam,mlod,mlom,melev, 1 ity,lerr) c fortran program to read inventory files on the prime c and build an interactive file database. c adn 4/29/87 durant, ok c read from SDATA c mst state number (missing first digit) c mstat station code c ndis c nam station name c mby beginning year of reacor c mbm beginning month of record c mey ending year of record c mem ending month of record c mlad latitude degrees c mlam latitude minutes c mlod longitude degrees c mlom longitude minutes c melev elevation (tens of feet) c ity type (1..4; used for maximum rainfall intensity) c county county common /blkx/idx(91,2) character*23 nam lerr=1 c deh print *,' getsta' 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 c*************** c JULIN c*************** 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 .eq. nyer) ii(2) = 29 do 30 j = 1,12 jday = jday -ii(j) if (jday.le.0) then nday = jday +ii(j) go to 35 end if 30 continue 35 jday = kday mo = j return end c********************************* c INTSR c********************************* subroutine intsr (tlat,tlon) c Routine intsr 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. c Opens and reads unit 13, file "statparm" c Calls lltoalbr c distance 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),dpx(12),srx(12),sdx(12),apx(12) 101 format (a32,7x,2f6.2,////) print *,' ************** intsr -- solar radiation' 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) c Find ten closest solar radiation and max .5 P stations (lines 3 to 7) j=0 3 j=j+1 read(13,101,end=7) statin,alat,alon !reads lat and lon in degrees call distance (tlat,tlon,alat,alon,dist) do 4 i=10,1,-1 if(dist.ge.far(i)) then k=i+1 go to 5 end if c Removed go to and fixed with if then by DEH c if(dist.lt.far(i)) go to 4 c k=i+1 c go to 5 4 continue k=1 5 do 6 i=10,k,-1 if (i.gt.1) then far(i)=far(i-1) sta(i)=sta(i-1) xlat(i)=xlat(i-1) xlon(i)=xlon(i-1) indx(i)=indx(i-1) 6 end if 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)) c finds the three closest stations that make an enclosing triangle 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).eq.tx) then ! fuzz slpx1=9999. else slpx1=(y(l)-ty)/(x(l)-tx) end if cptxt1=y(l)-slpx1*x(l) if(x(k).eq.x(i)) then slpki=9999. else slpki=(y(k)-y(i))/(x(k)-x(i)) end if cptki=y(k)-slpki*x(k) if(slpx1.eq.slpki) then sx1=9999.0 else sx1=(cptki-cptxt1)/(slpx1-slpki) end if sy1=slpx1*sx1+cptxt1 if((sx1.ge.x(k).and.sx1.le.x(i)).or. 1 (sx1.le.x(k).and.sx1.ge.x(i))) then if((sy1.ge.y(k).and.sy1.le.y(i)).or. 1 (sy1.le.y(k).and.sy1.ge.y(i))) then c potential problem here -- if x(k) == tx, then cptxt2=y(k)-9999.*x(k) if(x(k).eq.tx) then slpx2=9999. else slpx2=(y(k)-ty)/(x(k)-tx) end if cptxt2=y(k)-slpx2*x(k) if(x(i).eq.x(l)) then slpil=9999. else slpil=(y(i)-y(l))/(x(i)-x(l)) end if 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. 1 (sx2.le.x(i).and.sx2.ge.x(l))) then if((sy2.ge.y(i).and.sy2.le.y(l)).or. 1 (sy2.le.y(i).and.sy2.ge.y(l))) then ic(1)=l ic(2)=k ic(3)=i go to 22 endif endif endif endif 21 continue 22 continue c If no triangulation found, then ic(1),ic(2),ic(3) equal zero. c c Checking for stations triangulated with extreme distance. if (ic(3).gt.0) then if (far(ic(3)).gt. 2000.) then print *,' WARNING! closest selected stations are: ', 1 far(ic(1)),far(ic(2)),far(ic(3)) ic(1)=0 ! ic(1)=0 to make program skip and use two stations end if endif c Finding weighting for triangulation of three stations. if(ic(1).ne.0) then 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) else c Continue on finding two stations to use if no triangulation found. 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.eq.0.0) then slpnor=9999. else slpnor=-1./slp end if 30 antnor=ty-slpnor*tx ssx=(antnor-antrcp)/(slp-slpnor) ssy=slpnor*ssx+antnor ssxtx=ssx-tx ssyty=ssy-ty d=far(k)+far(l)+((ssxtx*ssxtx+ssyty*ssyty)**0.5)/500. if(d.lt.dis1) then dis1=d ic(1)=l ic(2)=k end if 26 continue c only two points being used (could not find third station to triangulate) l=ic(1) k=ic(2) i=0 ! deh added so random station not reported w=far(l)+far(k) wgt(1)=(w-far(l))/w wgt(2)=(w-far(k))/w wgt(3)=0.0 endif c APPEND c Write stations and weights to triang.tmp to append to par file. DLS open(22,file="triang.tmp", status="unknown") if (i.gt.0) then write(22,77) sta(l),wgt(1),sta(k),wgt(2),sta(i),wgt(3) else write(22,78) sta(l),wgt(1),sta(k),wgt(2) end if 77 format(a19,f6.3,2(2x,a19,f6.3)) 78 format(a19,f6.3,2x,a19,f6.3) 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) if (i.gt.0) 1 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 (i.eq.0) then if (j.ne.indx(l).and.j.ne.indx(k))then read(13,'(a5,////)',end=29) aa go to 51 else read(13,104,end=29) tp5x,tp6x 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) endif else if (j.ne.indx(l).and.j.ne.indx(k).and.j.ne.indx(i))then read(13,'(a5,////)',end=29) aa go to 51 else read(13,104,end=29) tp5x,tp6x 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) endif end if 104 format(59x,2f6.0) 105 format (8x,12f6.0) if(j.eq.indx(l)) awgt=wgt(1) if(j.eq.indx(k)) awgt=wgt(2) if(i.ne.0 .and. 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********************************* c INTDP c********************************* subroutine intdp (tlat,tlon) c Subroutine "intdp" is just like "intsr" except c this one was added to handle dewpoint from a different file. c Opens and reads unit 13, file 'dewpoint" c Calls lltoalbr c distance c Note: unused symbol APX c unused symbol SDX c unused symbol SRX 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),dpx(12) ! ,srx(12),sdx(12),apx(12) 101 format (a29,38x,2f3.0,f4.0,f3.0) print *,' ************** intdp: dewpoint' 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) c find ten closest dewpoint reference stations (lines 3 to 7) 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.ge.far(i)) then k=i+1 go to 5 end if c Removed go to and fixed with if then by DEH c if(dist.lt.far(i)) go to 4 c k=i+1 c go to 5 4 continue k=1 5 do 6 i=10,k,-1 if (i.gt.1) then far(i)=far(i-1) sta(i)=sta(i-1) xlat(i)=xlat(i-1) xlon(i)=xlon(i-1) indx(i)=indx(i-1) 6 end if far(k)=dist xlat(k)=alat xlon(k)=alon sta(k)=statin indx(k)=j go to 3 c Converts latitude and longitude of closest stations to Albers projection 7 do 8 i=1,10 8 call lltoalbr (xlat(i),xlon(i),x(i),y(i)) c Finds three of closest stations that make triangle by intersection of lines. 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).eq.tx) then slpx1=9999. else slpx1=(y(l)-ty)/(x(l)-tx) endif cptxt1=y(l)-slpx1*x(l) if(x(k).eq.x(i)) then slpki=9999. else slpki=(y(k)-y(i))/(x(k)-x(i)) endif cptki=y(k)-slpki*x(k) if(slpx1.eq.slpki) then sx1=9999.0 else sx1=(cptki-cptxt1)/(slpx1-slpki) endif sy1=slpx1*sx1+cptxt1 if((sx1.ge.x(k).and.sx1.le.x(i)).or. 1 (sx1.le.x(k).and.sx1.ge.x(i))) then if((sy1.ge.y(k).and.sy1.le.y(i)).or. 1 (sy1.le.y(k).and.sy1.ge.y(i))) then if(x(k).eq.tx) then slpx2=9999. else slpx2=(y(k)-ty)/(x(k)-tx) endif cptxt2=y(k)-slpx2*x(k) if(x(i).eq.x(l)) then slpil=9999. else slpil=(y(i)-y(l))/(x(i)-x(l)) endif 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. 1 (sx2.le.x(i).and.sx2.ge.x(l))) then if((sy2.ge.y(i).and.sy2.le.y(l)).or. 1 (sy2.le.y(i).and.sy2.ge.y(l))) then ic(1)=l ic(2)=k ic(3)=i go to 22 endif endif endif endif 21 continue 22 continue C If no triangulation found, then ic(1),ic(2),ic(3) equal zero. c c Checking for stations triangulated with extreme distance. if (ic(3).gt.0 .and. far(ic(3)).gt. 2000.) then print *,' WARNING! closest selected stations are: ', 1 far(ic(1)),far(ic(2)),far(ic(3)) ic(1)=0 end if if(ic(1).ne.0) then 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) c deh c if (wgt(3).lt.0.01) then c ic(3)=1 c i=ic(3) c w=far(l)+far(k) c wgt(1)=(w-far(l))/w c wgt(2)=(w-far(k))/w c wgt(3)=wgt(1) c end if else dis1=99999999. do 26 k=2,10 do 26 l=1,k-1 if(x(l).eq.x(k)) then slp=9999. else slp=(y(l)-y(k))/(x(l)-x(k)) endif antrcp=y(l)-slp*x(l) if(slp.eq.0.0) then slpnor=9999. else slpnor=-1./slp endif antnor=ty-slpnor*tx ssx=(antnor-antrcp)/(slp-slpnor) ssy=slpnor*ssx+antnor ssxtx=ssx-tx ssyty=ssy-ty d=far(k)+far(l)+((ssxtx*ssxtx+ssyty*ssyty)**0.5)/500. if(d.lt.dis1) then dis1=d ic(1)=l ic(2)=k endif 26 continue l=ic(1) k=ic(2) c deh "i" was "random" if third weight is zero; c station on screen didn't match what is in file. c par file had station 3 = station 1, so we do the same here. c ic(3)=1 ic(3) not used in subroutine beyond this point DLS c i=ic(3) i=0 w=far(l)+far(k) wgt(1)=(w-far(l))/w wgt(2)=(w-far(k))/w wgt(3)=0.0 c wgt(3)=wgt(1) end if c Write stations and weights to triang.tmp to append to par file. DLS if (i.gt.0) then write(22,77) sta(l),wgt(1),sta(k),wgt(2),sta(i),wgt(3) else write(22,78) sta(l),wgt(1),sta(k),wgt(2) end if 77 format(a19,f6.3,2(2x,a19,f6.3)) 78 format(a19,f6.3,2x,a19,f6.3) 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) if (i.gt.0) 1 write(*,102) indx(i),sta(i),xlat(i),xlon(i),far(i),wgt(3) print *,' The 10 closest stations are: ',(sta(ii),ii=1,10) 102 format(i4,a22,4f9.3) close(13) open (13,file='dewpoint',status='old') j=0 51 j=j+1 if (i.eq.0) then if (j.ne.indx(l).and.j.ne.indx(k)) then read(13,'(a5)',end=29) aa go to 51 else read(13,'(31x,12f3.0)',end=29) (dpx(m),m=1,12) endif else if (j.ne.indx(l).and.j.ne.indx(k).and.j.ne.indx(i)) then read(13,'(a5)',end=29) aa go to 51 else read(13,'(31x,12f3.0)',end=29) (dpx(m),m=1,12) endif end if if(j.eq.indx(l)) awgt=wgt(1) if(j.eq.indx(k)) awgt=wgt(2) if(i.ne.0 .and. 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********************************* c INTTP c********************************* subroutine inttp (tlat,tlon,tp) c Subroutine inttp was written by Gene Gander to interpolate values at c specified Lat-Lon points using closest enclosing triangles of three c points where there exist known values for the parameter of interest. c In the case of this subroutine, the parameter of interest is the c probability curve of time to peak rain intensity. c Opens and reads unit 13, file "timepeak.tot" c Calls lltoalbr c distance 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),tp(12),tpx(12) c 101 format (6x,a32,1x,4f5.0) 101 format (13x,a25,1x,4f5.0) print *,' ************** inttp -- time to peak' do 61 i=1,12 tp(i)=0.0 61 tpx(i)=0.0 1 do 2 i=1,10 2 far(i)=99999999.0 open(13,file='timepeak.tot',status='old') call lltoalbr (tlat,tlon,tx,ty) c Find ten closest time to peak stations j=0 3 j=j+1 read(13,101,end=7) statin,alatd,alatm,alond,alonm read(13,'(a)') aa 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.ge.far(i)) then k=i+1 go to 5 end if 4 continue k=1 5 do 6 i=10,k,-1 if (i.gt.1) then far(i)=far(i-1) sta(i)=sta(i-1) xlat(i)=xlat(i-1) xlon(i)=xlon(i-1) indx(i)=indx(i-1) end if 6 continue 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).eq.tx) then slpx1=9999. else slpx1=(y(l)-ty)/(x(l)-tx) end if cptxt1=y(l)-slpx1*x(l) if(x(k).eq.x(i)) then slpki=9999. else slpki=(y(k)-y(i))/(x(k)-x(i)) end if cptki=y(k)-slpki*x(k) if(slpx1.eq.slpki) then sx1=9999.0 else sx1=(cptki-cptxt1)/(slpx1-slpki) end if sy1=slpx1*sx1+cptxt1 if((sx1.ge.x(k).and.sx1.le.x(i)).or. 1 (sx1.le.x(k).and.sx1.ge.x(i))) then if((sy1.ge.y(k).and.sy1.le.y(i)).or. 1 (sy1.le.y(k).and.sy1.ge.y(i))) then if(x(k).eq.tx) then slpx2=9999. else slpx2=(y(k)-ty)/(x(k)-tx) end if cptxt2=y(k)-slpx2*x(k) if(x(i).eq.x(l)) then slpil=9999. else slpil=(y(i)-y(l))/(x(i)-x(l)) end if 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. 1 (sx2.le.x(i).and.sx2.ge.x(l))) then if((sy2.ge.y(i).and.sy2.le.y(l)).or. 1 (sy2.le.y(i).and.sy2.ge.y(l))) then ic(1)=l ic(2)=k ic(3)=i go to 22 endif endif endif endif 21 continue 22 continue if (ic(3).gt.0 .and. far(ic(3)).gt. 2000) then print *,' WARNING! closest selected stations are: ', 1 far(ic(1)),far(ic(2)),far(ic(3)) ic(1)=0 endif if (ic(1).ne.0) then 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) else dis1=99999999. do 26 k=2,10 do 26 l=1,k-1 if(x(l).eq.x(k)) then slp=9999. else slp=(y(l)-y(k))/(x(l)-x(k)) endif antrcp=y(l)-slp*x(l) if(slp.eq.0.) then slpnor=9999. else slpnor=-1./slp end if antnor=ty-slpnor*tx ssx=(antnor-antrcp)/(slp-slpnor) ssy=slpnor*ssx+antnor ssxtx = ssx-tx ssyty = ssy-ty d=far(k)+far(l)+((ssxtx*ssxtx+ssyty*ssyty)**0.5)/500. if(d.lt.dis1) then dis1=d ic(1)=l ic(2)=k endif 26 continue l=ic(1) k=ic(2) c deh "i" was "random" if third weight is zero; c station on screen didn't match what is in file. c File had station 3 = station 1, so we do the same here. c ic(3)=1 !ic(3) not used in subroutine beyond this point DLS c i=ic(3) i=0 w=far(l)+far(k) wgt(1)=(w-far(l))/w wgt(2)=(w-far(k))/w wgt(3)=0.0 c wgt(3)=wgt(1) ! **************************** end if c Write stations and weights to triang.tmp to append to par file. DLS if (i.gt.0) then write(22,77) sta(l),wgt(1),sta(k),wgt(2),sta(i),wgt(3) else write(22,78) sta(l),wgt(1),sta(k),wgt(2) end if 77 format(a19,f6.3,2(2x,a19,f6.3)) 78 format(a19,f6.3,2x,a19,f6.3) close(22) 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) if (i.ne.0) 1 write(*,102) indx(i),sta(i),xlat(i),xlon(i),far(i),wgt(3) 102 format(i4,2x,a32,4f9.3) close(13) open (13,file='timepeak.tot',status='old') c print *,'looking for records',indx(l) c print *,' ',indx(k) c if (i.gt.0) print *,' ', indx(i) j=0 51 j=j+1 read(13,'(a5)',end=29) aa if (i.eq.0) then if (j.ne.indx(l).and.j.ne.indx(k)) then read(13,'(a5)') aa go to 51 else read(13,'(7x,12f6.3)',end=29) (tpx(m),m=1,12) endif else if (j.ne.indx(l).and.j.ne.indx(k).and.j.ne.indx(i)) then read(13,'(a5)') aa go to 51 else read(13,'(7x,12f6.3)',end=29) (tpx(m),m=1,12) endif end if if (j.eq.indx(l)) awgt=wgt(1) if (j.eq.indx(k)) awgt=wgt(2) if (i.ne.0 .and. j.eq.indx(i)) awgt=wgt(3) do 52 n=1,12 tp(n)= tp(n)+awgt*tpx(n) c print *,n,tp(n),awgt 52 continue go to 51 29 close(13) return end c********************************* c INTWIND c********************************* subroutine intwind (tlat,tlon,wnd,wgt,site) c Opens and reads unit 13, file "idxall" c Opens and reads unit 14, three nearest stations, file sta(iii) c Note: unused symbol AA c character*29 sta(15),statin ! depends on path length in data file! character*16 sta(15),statin c character*5 aa character*19 site(3) dimension wnd(12,4,17),wgt(3) dimension xlat(15),xlon(15),x(15),y(15),far(15),indx(15),ic(3) 1 ,temp(12,4,17) c deh 101 format (a29,29x,2f3.0,f4.0,f3.0) ! ditto 101 format (a16,29x,2f3.0,f4.0,f3.0) do 61 i=1,12 do 61 j=1,4 do 61 k=1,17 61 wnd(i,j,k)=0.0 1 do 2 i=1,10 2 far(i)=99999999.0 c deh open(13,file='/usr3/c/wind/idxall',status='old') open(13,file='idxall',status='old') print *,' ************** intwind -- wind' call lltoalbr (tlat,tlon,tx,ty) c Find ten closest wind stations 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.ge.far(i)) then k=i+1 go to 5 end if 4 continue k=1 5 do 6 i=10,k,-1 if (i.gt.1) then far(i)=far(i-1) sta(i)=sta(i-1) xlat(i)=xlat(i-1) xlon(i)=xlon(i-1) indx(i)=indx(i-1) end if 6 continue 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).eq.tx) then slpx1=9999. else slpx1=(y(l)-ty)/(x(l)-tx) end if cptxt1=y(l)-slpx1*x(l) if(x(k).eq.x(i)) then slpki=9999. else slpki=(y(k)-y(i))/(x(k)-x(i)) end if cptki=y(k)-slpki*x(k) if(slpx1.eq.slpki) then sx1=9999.0 else sx1=(cptki-cptxt1)/(slpx1-slpki) end if sy1=slpx1*sx1+cptxt1 if((sx1.ge.x(k).and.sx1.le.x(i)).or. 1 (sx1.le.x(k).and.sx1.ge.x(i))) then if((sy1.ge.y(k).and.sy1.le.y(i)).or. 1 (sy1.le.y(k).and.sy1.ge.y(i))) then if(x(k).eq.tx) then slpx2=9999. else slpx2=(y(k)-ty)/(x(k)-tx) end if cptxt2=y(k)-slpx2*x(k) if(x(i).eq.x(l)) then slpil=9999. else slpil=(y(i)-y(l))/(x(i)-x(l)) end if 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. 1 (sx2.le.x(i).and.sx2.ge.x(l))) then if((sy2.ge.y(i).and.sy2.le.y(l)).or. 1 (sy2.le.y(i).and.sy2.ge.y(l))) then ic(1)=l ic(2)=k ic(3)=i go to 22 endif endif endif endif 21 continue 22 continue if (ic(3).gt.0 .and. far(ic(3)).gt. 2000) then print *,' WARNING! closest selected stations are: ', 1 far(ic(1)),far(ic(2)),far(ic(3)) ic(1)=0 site(3)=' ' endif if(ic(1).ne.0) then l=ic(1) k=ic(2) i=ic(3) if(far(k).lt.0.001.and.far(l).lt.0.001) then wgt(1)=0.5 wgt(2)=0.5 wgt(3)=0.0 else 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) end if else dis1=99999999. do 26 k=2,10 do 26 l=1,k-1 if(x(l).eq.x(k)) then slp=9999. else slp=(y(l)-y(k))/(x(l)-x(k)) endif antrcp=y(l)-slp*x(l) if(slp.eq.0.) then slpnor=9999. else slpnor=-1./slp end if antnor=ty-slpnor*tx ssx=(antnor-antrcp)/(slp-slpnor) ssy=slpnor*ssx+antnor ssxtx = ssx-tx ssyty = ssy-ty d=far(k)+far(l)+((ssxtx*ssxtx+ssyty*ssyty)**0.5)/500. if(d.lt.dis1) then dis1=d ic(1)=l ic(2)=k endif 26 continue l=ic(1) k=ic(2) ic(3)=0 c deh "i" was "random" if third weight is zero; c station on screen didn't match what is in file. c File had station 3 = station 1, so we do the same here. c i=ic(3) c deh i=0 w=far(l)+far(k) wgt(1)=(w-far(l))/w wgt(2)=(w-far(k))/w wgt(3)=0.0 site(3)=' ' c print *, 'wgt(3)=0' c print *, 'l,k,i=', l,k,i c print *, 'ic(1,2,3)=',ic(1),ic(2),ic(3) end if 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) if(i.gt.0) 1 write(*,102) indx(i),sta(i),xlat(i),xlon(i),far(i),wgt(3) print *,' The closest stations are: ' print *, sta(1), far(1) print *, sta(2), far(2) print *, sta(3), far(3) 102 format(i4,a30,4f9.3) close(13) do 29 mm=1,3 iii=ic(mm) if (iii.ne.0) then open (14,file=sta(iii),status='old',err=2222) read(14,'(a)') site(mm) do 71 m1=1,16 do 71 m2=1,4 71 read(14,108) (temp(nn,m2,m1),nn=1,12) read(14,108) (temp(nn,1,17),nn=1,12) 108 format(12f6.2) do 52 n1=1,12 do 52 n2=1,4 do 52 n3=1,17 52 wnd(n1,n2,n3)= wnd(n1,n2,n3)+wgt(mm)*temp(n1,n2,n3) close(14) end if 29 continue goto 2223 2222 write(*,*) ' Error opening ',sta(iii) stop 2223 continue return end c********************************* c DISTANCE c********************************* subroutine distance (alat,alon,blat,blon,dist) c Routine distance 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. 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. 1 abs(alat).ge.abs(blat-.01).and. 2 abs(alon).le.abs(blon+.01).and. 3 abs(alon).ge.abs(blon-.01)) then dist=0.0 else argue=cos(a)*cos(c)*cos(d-b)+sin(a)*sin(c) if(argue.gt.1.0) argue=1.0 if(argue.lt.-1.0) argue=-1.0 dist=radius*acos(argue) endif return end c********************************* c LLTOALBR c********************************* SUBROUTINE LLTOALBR(alat,alon,x,y) 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. c Standard parallels of 29.5 deg and 45.5 deg seem to be quite c standard for continental US work. [gis.state.ga.us; www.ce.utexas.edu] c (although ftp.sas.com uses 0 and 45) c For South Africa, 18S and 32S are used c [http://www.ccwr.ac.za/~/lynch2/coords.html] c The consensus sems to be 63710022.77 meters for Earth radius c [www.ccwr.ac.za; c "The Albers Equal Area projection has the property that the area c bounded by any pair of parallels and meridians is exactly reproduced c between the image of those parallels and meridians in the projected c domain, that is, the projection preserves the correct area of the c earth though distorts direction, distance, and shape somewhat." c http://www.ce.utexas.edu/prof/maidment/CE397/ex397/ex397.htm c Exercise 3: Map Projections [7/26/1999] c Does it make sense to use the Albers Equal-Area projection when we c are interested in direction and distance rather than area? (Also, c it is limited to the Continental US unless parameters are set for c other areas). Look into the Robinson projection, "...a relatively c new style of map projection for the earth designed to present the c whole earth witha minimum of distortion at any location." [ibid] c The coordinates are in meters in the projected coordinate system. c 29 deg 30 min (29.5 deg) and 45 deg 30 min (45.5 deg) -96 23 c produce a reasonable map for the continental United States, but c are not suitable for Alaska and Hawaii (They are rotated maybe c 75 degrees clockwise). c For regions that are primarily East-West in extent, the Albers c Equal Area or Lambert Conformal Conic are good choices for map c projections. The Central Meridian should be through the middle of c the region of interest, the Reference Latitude should be wherever c you think the center of a coordinate system should be (usually at c the center or below the bottom of the extent of the geographic c features. The two standard parallels should be located approximately c 1/6 from the top of the geographic extent of the mapped features. c "Albers Projection is a conic projection which uses two standard c parallels to reduce some of the distortion produced when only one c standard parallel is used. Although neither shape nor linear scale c are truly correct, the distortion of these properties is minimized c in the region between the standard parallels. c http://gis.state.ga.us/gisfaq/proj.html c Projections and Datums FAQ [7/26/1999] c (Assuming the NAD83 datum) for Standard Parallel 1 = 29.5 deg N, c Standard Parallel 2 = 45.5 deg N, a Central Meridian of 97 degrees c West (-97) a Reference Latitude of 23 degrees North, and zero c Xshift and Yshift, a point at (97W,23N) degrees will become (0,0) c meters in the projected Albers coordinate system. [www.ce.utexas.edu] REAL lamda0,lambda pi=3.141592654 sp1 = 29.5 ! standard parallel 1 (29 degrees 30 min N) sp2 = 45.5 ! standard parallel 2 (45 degrees 30 min N) cmer = -97.0 ! central meridian (97 degrees West) r=6370997. ! earth radius, in meters phi0=90.*PI/180. phi1=sp1*PI/180. phi2=sp2*PI/180. lamda0=cmer*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 C PROGRAM TO CONVERT COORDINATES IN DEGREES TO COORDINATES IN METRES C BY USING THE ALBERS EQUAL AREA TRANSFORMATION. c from ftp://water.ccwr.ac.za/pub/lynch/coords/albers.for [7/26/1999] C c Only the returned arguments, i.e. answers in metres must be double precision c in the calling program. c c SUBROUTINE ALBERS(XLAT,XLONG,ZYP,ZXP) c IMPLICIT REAL*8 (Z) c Enter LAT (+ LAT) and LONG in minutes of a degree ' c For example c LAT = 1800 minutes (30 degrees) c LONG = 1500 minutes (25 degrees) c answer c our X = 95929.280 metres c our Y = -3281845.424 metres south of the equator c ZLAT =DBLE(ABS(XLAT)) c ZLONG=DBLE(ABS(XLONG)) c ZLAT = ZLAT*DACOS(-1.0D0)/180.0D0/60.0D0 c ZLONG= ZLONG*DACOS(-1.0D0)/180.0D0/60.0D0 c ZP2 = 32.0D0*DACOS(-1.0D0)/180.0D0 c ZP1 = 18.0D0*DACOS(-1.0D0)/180.0D0 c ZXL0= 24.0D0*DACOS(-1.0D0)/180.0D0 c ZYL0= 0.0D0 c ZSCALE=1000.0D0 c ZR = 6371002.77D0/ZSCALE/1000.0D0 c ZAK = (DSIN(ZP1) + DSIN(ZP2))/2.0D0 c ZC = DCOS(ZP1)*DCOS(ZP1) + 2.0D0*ZAK*DSIN(ZP1) c ZR0 = ZR*DSQRT(ZC-2.0D0*ZAK*DSIN(ZYL0))/ZAK c ZRP = ZR*DSQRT(ZC-2.0D0*ZAK*DSIN(ZLAT))/ZAK c ZPP = ZAK*(ZLONG-ZXL0) c ZXP = -(ZR0-ZRP*DCOS(ZPP))*1000.0D0*ZSCALE c ZYP = ZRP*DSIN(ZPP)*1000.0D0*ZSCALE c RETURN c END END