c ******************************************************************************* c * ****** program dat2par2 ***** c c Create CLIGEN input parameter ".par" file from various input files: c c *.dat precipitation and temperature data for the station c stations.dat list of station numbers, names, geolocations c statparm solar radiation, MAX .5 P, and dew point data for triangulation c dewpoint dewpoint data for triangulation c timepeak.tot time to peak distributions data for triangulation c idxall index of wind files with latitude and longitude c xx/*.par one directory per state of wind .par files c c Modification of Calparms.f77 found on computer of c Arlin Nicks c U.S.D.A. Agricultural Research Service c by c David Hall (DEH) and Dayna Scheele (DLS) c U.S.D.A. Forest Service, Moscow Forestry Sciences Lab c 06/1999 DEH, DLS c display non-random third-closest station if only 2 are used. c remove computed gotos from subroutines. c Documentation added, code cleaned up a bit 08/2000 DEH c remove extraneous commons and out-of-date comment lines c c ****************************************************************** c c Main opens, reads, closes unit 16 'list' c opens, reads, closes unit 7 '*.dat' (fil) c opens, reads, closes unit 11 'stations.dat' c opens, reads, closes unit 22 'triang.tmp' c opens, writes, closes unit 12 '*.par' (bfil) c Intsr opens, reads, closes unit 13 'statparm' c opens, writes unit 22 'triang.tmp' c Intdp opens, reads, closes unit 13 'dewpoint' c writes unit 22 'triang.tmp' c Inttp opens, reads, closes unit 13 'timepeak.tot' c writes, closes unit 22 'triang.tmp' c Intwind opens, reads, closes unit 13 'idxall' c opens, reads, closes unit 14 '*.par' c Calls julin -- given day of year return month and day of month 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 distance -- returns distance between two points of latitude-longitude c lltoalbr -- convert lat-lon points to X-Y coordinates !MS$DEBUG common /blksr/ solrad(12),sdsol(12),ax5p(12),dewpt(12),tp5,tp6 dimension 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) dimension sumx(12),ssx(12),sd(12),xmean(12),sk(12),sx3(12) dimension pmax(12),title(80),nc(12),misda(12),nrec(12) dimension wc(12),rc(12),dc(12),tc(12),pww(12),pdw(12),pdd(12), 1 pwd(12),xmonth(12),ngrec(12),ntrec(12) C dimension rst(12,3),prw(2,12),k1(4),k2(4),mo,v1,v2,st,lw,nrs dimension rst(12,3),prw(2,12),k1(4),k2(4) dimension sumt(2,12),sst(2,12),sdt(2,12),xmeant(2,12), 1 nrtc(2,12),misdt(2,12) dimension nr(366),tp(12),wnd(12,4,17),wgt(3) ! ,prc(366) integer mo real*8 st character*32 fil,bfil ! ,afil character*12 stnum character*23 nam character*1 b(32) c character*15 county character*80 line character*19 site(3) character*4 el(3),elem character*7 point(17),stat(3) 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 open(16,file='list',status='old',access='sequential',err=2050) c == get name of .dat file to process == 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 c == build output file name == 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 if (stnum(1:1) .eq. '#') goto 450 c first line after # is skipped 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 == found match in stations.dat == c == open .dat file == open (7,file=fil,status='old',access='sequential',err=454) goto 455 454 print *,' Data file not found for station' goto 10 c == create output .par file == 455 open(12,file=bfil,status='unknown') do 40 m=1,12 sumx(m)=0.0 ssx(m)=0.0 sx3(m)=0.0 pmax(m)=0.0 xmean(m)=0.0 xmonth(m)=0.0 nwc(m)=0 ndc(m)=0 nrc(m)=0 nrec(m)=0 ntrec(m)=0 ntc(m)=0 misda(m)=0 do 40 j=1,2 sumt(j,m)=0.0 sst(j,m)=0.0 sdt(j,m)=0.0 xmeant(j,m)=0.0 misdt(j,m)=0 nrtc(j,m)=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 if (elem.ne.el(1)) go to 210 ! look for 'PRCP' record l=1 years=years+1.0 ! flip the year counter lda = 365 if (((ny/4*4)-ny).eq.0) lda = 366 ! Add a day for nominal leap year c == for each day in the record (one record per year) do 200 i = 1,lda call julin(mo,ida,ny,i) if (nr(i).eq.9999) misda(mo)=misda(mo)+1 ! count missing days ntc(mo)=ntc(mo)+1 ! count days for month if (nr(i).gt.0.and.nr(i).lt.9999) nrc(mo)=nrc(mo)+1 ! # non-missing days if (nr(i).le.0.and.nt.ne.0) ndc(mo)=ndc(mo)+1 ! # wet -> dry if (nr(i).gt.0.and.nr(i).ne.9999.and.nt.gt.0)nwc(mo)=nwc(mo)+1 ! # wet -> wet nt=nr(i) ! nt = nr(i-1) if (nr(i).eq.9999) nt=0 ! nt = 0 (dry) if missing data if (nr(i).ne.0 .and. nr(i).lt.9998) then st=nr(i) ! precip in hundredths of inches st=st*.01 ! precip in inches x(mo,ida)=st sumx(mo)=sumx(mo)+st ! sum of precip for wet days for month ssx(mo)=ssx(mo)+(st*st) ! sum of precip**2 for wet days for month sx3(mo)=sx3(mo)+(st*st*st) ! sum of precip**3 for wet days if (st.gt.pmax(mo)) pmax(mo)=st ! gather maximum precip value end if 200 continue do 201 m=1,12 201 nrec(m)=nrec(m)+1 go to 1 ! get another record (year) of data 210 ld=1 c lda=365 c if (((ny/4*4)-ny).eq.0) lda=366 j=1 ! 'TMAX' record if (elem.eq.el(3)) j=2 ! 'TMIN' record do 220 i=1,lda call julin(mo,ida,ny,i) if (nr(i).eq.9999) misdt(j,mo)=misdt(j,mo) + 1 ! count missing days c if (nr(i).lt.9999) nrtc(j,mo)=nrtc(j,mo) + 1 ! count non-missing days if (nr(i).ne.9999) then nrtc(j,mo)=nrtc(j,mo) + 1 ! count non-missing days st=nr(i) ! temperature sumt(j,mo)=sumt(j,mo)+st ! sum of temperatures for month sst(j,mo)=sst(j,mo)+(st*st) ! sum of temperatues**2 for month end if 220 continue go to 1 ! try for another record (year) c == came to end of .dat file == c == calculate p(w|w), p(d|w), p(w|d), p(d|d) == 50 do 250 m=1,12 wc(m)=nwc(m) rc(m)=nrc(m) dc(m)=ndc(m) tc(m)=ntc(m) pww(m)=wc(m)/rc(m) pdw(m)=1.-pww(m) pwd(m)=(rc(m)-wc(m))/(tc(m)-rc(m)) pdd(m)=1.-pwd(m) ! was: pdd(m)=1.-pdw(m). [= pww(m)]. Not used. 250 continue c == calculate std dev and skew coefficients for precip == do 300 m=1,12 if (nrc(m).lt.3) nrc(m)=3 an = nrc(m) ! number non-missing for month am = nrec(m) ! number of entries (days) for month xmonth(m)=sumx(m)/am ! avg precip per day for month xmean(m)=sumx(m)/an ! avg precip per day for non-missing days sd(m)=ssx(m)-(sumx(m)*sumx(m)/an) sd(m)=sqrt(sd(m)/(an-1.0)) sum3=sumx(m)*sumx(m)*sumx(m) sk(m)=an*an*sx3(m)-3.0*an*sumx(m)*ssx(m)+2.0*sum3 sk(m)=sk(m)/(an*(an-1.0)*(an-2.0)) 300 sk(m)=sk(m)/(sd(m)*sd(m)*sd(m)) c == calculate std dev for temperatures == do 310 j=1,2 do 310 m=1,12 an=nrtc(j,m) xmeant(j,m)=sumt(j,m)/an sdt(j,m)=sst(j,m)-(sumt(j,m)*sumt(j,m)/an) sdt(j,m)=sqrt(sdt(j,m)/(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)(m,nrec(m),misda(m),nrc(m),nwc(m),ndc(m),ntc(m), 1 xmonth(m),pmax(m),m=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)(m,pww(m),pwd(m),xmean(m),sd(m),sk(m),m=1,12) 103 format(1x,i4,f10.3,f8.3,f7.3,f8.3,f8.3) do 130 m=1,12 xmonth(m)=0.0 sumx(m)=0.0 ngrec(m)=0 prw(1,m)=pww(m) prw(2,m)=pwd(m) rst(m,1)=xmean(m) rst(m,2)=sd(m) 130 rst(m,3)=sk(m) write(*,1030) (xmeant(1,m),m=1,12) write(*,1030) (xmeant(2,m),m=1,12) write(*,1030) (sdt(1,m),m=1,12) write(*,1030) (sdt(2,m),m=1,12) write(*,1031) (misdt(1,m),m=1,12) write(*,1031) (misdt(2,m),m=1,12) write(*,1031) (nrtc(1,m),m=1,12) write(*,1031) (nrtc(2,m),m=1,12) 1030 format(1x,12f6.2) 1031 format(1x,12i6) ity = itype print *, nst,' ',nstat,' ',nam 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(m),m=1,12) 4030 format(' MEAN P ',12f6.2) write(12,4040)(sd(m),m=1,12) 4040 format(' S DEV P',12f6.2) write(12,4050)(sk(m),m=1,12) 4050 format(' SKEW P',12f6.2) write(12,4060)(pww(m),m=1,12) 4060 format(' P(W/W) ',12f6.2) write(12,4070)(pwd(m),m=1,12) 4070 format(' P(W/D) ',12f6.2) write(12,4080)(xmeant(1,m),m=1,12) 4080 format(' TMAX AV',12f6.2) write(12,4085)(xmeant(2,m),m=1,12) 4085 format(' TMIN AV',12f6.2) write(12,4090)(sdt(1,m),m=1,12) 4090 format(' SD TMAX',12f6.2) write(12,4095)(sdt(2,m),m=1,12) 4095 format(' SD TMIN',12f6.2) write(12,4100)(solrad(m),m=1,12) 4100 format(' SOL.RAD',12F6.0) write(12,4110)(sdsol(m),m=1,12) 4110 format(' SD SOL ',12F6.1) write(12,4120)(ax5p(m),m=1,12) 4120 format(' MX .5 P',12F6.2) c********************************************************************* c DEW POINT c********************************************************************* call intdp(xlat,xlon, dewpt) write(12,4130)(dewpt(m),m=1,12) 4130 format(' DEW PT ',12F6.2) c********************************************************************* c TIME TO PEAK c********************************************************************* call inttp(xlat,xlon,tp) write(12,4131)(tp(m),m=1,12) 4131 format('Time Pk ',12f6.3) c********************************************************************** c WIND DATA c********************************************************************** call intwind (xlat,xlon,wnd,wgt,site) do 4134 i=1,16 ! compass points write(12,4102) point(i),(wnd(m,1,i),m=1,12) do 4134 k=2,4 4134 write(12,4102) stat(k-1),(wnd(m,k,i),m=1,12) write(12,4102) point(17),(wnd(m,1,17),m=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 ! go get another .dat file 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 (12) stop end c*************** c JULIN c*************** subroutine julin (mo,nday,nyer,jday) c given Julian day (1..366) (jday) and nominal year (nyer) c return month (1..12) (mo) and day of month (1..31) (nday) c Leap year is encountered every 4 years only 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 m = 1,12 jday = jday -ii(m) if (jday.le.0) then nday = jday +ii(m) go to 35 end if 30 continue 35 jday = kday mo = m return end c********************************* c INTSR c********************************* subroutine intsr (tlat,tlon) c Interpolate values at specified lat-lon points using closest enclosing c triangles of three points where there exist known values for the c parameters of interest c (monthly mean and standard deviation of solar radiation) c Written Nov.,1988 by Gene Gander c Modified 1999 DLS and DEH c Intsr opens, reads, closes unit 13 'statparm' c opens, writes unit 22 'triang.tmp' 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 m=1,12 solrad(m)=0.0 sdsol(m)=0.0 ax5p(m)=0.0 dewpt(m)=0.0 srx(m)=0.0 sdx(m)=0.0 apx(m)=0.0 61 dpx(m)=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) 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)) 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 m=1,12 solrad(m)= solrad(m)+awgt*srx(m) sdsol(m)= sdsol(m)+awgt*sdx(m) ax5p(m)= ax5p(m)+awgt*apx(m) 52 dewpt(m)= dewpt(m)+awgt*dpx(m) tp5=tp5+awgt*tp5x tp6=tp6+awgt*tp6x go to 51 29 close(13) return end c********************************* c INTDP c********************************* subroutine intdp (tlat,tlon,dewpt) c Interpolates values at specified Lat-Lon points using closest enclosing c triangles of three points where there exist known monthly values for dewpoint c Gene Gander c Modified DLS and DEH 2000 c Intdp opens, reads, closes unit 13 'dewpoint' c writes unit 22 'triang.tmp' c Calls lltoalbr c distance c Note: unused symbol APX c unused symbol SDX c unused symbol SRX c common /blksr/ solrad(12),sdsol(12),ax5p(12),dewpt(12),tp5,tp6 dimension dewpt(12) 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 m=1,12 dewpt(m)=0.0 61 dpx(m)=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) end if 6 continue 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 m=1,12 52 dewpt(m)= dewpt(m)+awgt*dpx(m) go to 51 29 close(13) return end c********************************* c INTTP c********************************* subroutine inttp (tlat,tlon,tp) c Interpolates values at specified Lat-Lon points using closest enclosing c triangles of three points where there exist known values for the c probability curve of time to peak rain intensity. c --Gene Gander c Modified 2000 DLS and DEH c Inttp opens, reads, closes unit 13 'timepeak.tot' c writes, closes unit 22 'triang.tmp' 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 m=1,12 tp(m)=0.0 61 tpx(m)=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 m=1,12 tp(m)= tp(m)+awgt*tpx(m) c print *,n,tp(m),awgt 52 continue go to 51 29 close(13) return end c********************************* c INTWIND c********************************* subroutine intwind (tlat,tlon,wnd,wgt,site) c Interpolates values at specified Lat-Lon points using closest enclosing c triangles of three points where there exist known values for the c wind direction percantages and speeds c --Gene Gander c Modified 2000 DLS and DEH c Intwind opens, reads, closes unit 13 'idxall' c opens, reads, closes unit 14 '*.par' c 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 m=1,12 do 61 j=1,4 do 61 k=1,17 61 wnd(m,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 ! wind station number iii=ic(mm) if (iii.ne.0) then open (14,file=sta(iii),status='old',err=2222) ! wind station .par file read(14,'(a)') site(mm) do 71 m1=1,16 ! compass point do 71 m2=1,4 ! %, mean, SD, skew 71 read(14,108) (temp(m,m2,m1),m=1,12) ! month read(14,108) (temp(m,1,17),m=1,12) ! calm % row 108 format(12f6.2) do 52 m=1,12 ! month do 52 n2=1,4 ! parameter (%, mean...) do 52 n3=1,17 ! compass point 52 wnd(m,n2,n3)= wnd(m,n2,n3)+wgt(mm)*temp(m,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 Returns the distance (dist) between two points of latitude-longitude (degrees) c in the c same units as the earth's radius (presently equal to 6370.997 kilometers) c Written Nov. 1988, Gene Gander. c Modified 2000, DEH 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 CONVERT 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. c THE SUBROUTINE ASSUMES AN EARTH RADIUS OF 6370997 METERS AND A VERTICAL c MERIDIAN OF LONGITUDE OF 97 DEGREES. c IF THESE ASSUMPTIONS ARE NOT VALID THEN MAKE APPROPRIATE CHANGES. c Following comments gathered by DEH 2000 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 END