program geofour Implicit double precision (a-h,o-z) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c g e o f o u r c c program for gravity field modelling by fft. a input data grid c e.g. generated by 'geogrid' is transformed, using kernels given c in the frequency domain. in the simple mode (<10) a gravity or geoid c grid is transformed assuming data to be given on a common level. c if data is given on the surface of the topography a first-order c harmonic continuation scheme may be used, and in this case a c height grid must also be input. c c c input: c ------ c c c (or possibly dummy name) c c mode, attkm, c (hkm - for mode 6 only) c rfi1, rla1, inn, ine, iwndow c c c where c c is the grid file containing data grid in standard c (unit 20) format, i.e. scanned in e-w lines from n to s, c initiated by label (lat1,lat2,lon1,lon2,dlat,dlon) c describing grid limits. utm grids may be used, in c this case lat/lon should be northing/easting in c meter, followed by ellipsoid number (1:wgs84, c 2:ed50, 3:nad27) and utm zone. c a subgrid of the data grid is analyzed, see below. c c is a height grid file on text format corresponding c (unit 21) to . not needed in the simple mode (dummy name c must be specified). the height grid must have the c same spacings and relative position as the data c grid, and must have at least the wanted area in common. c c outputfile (output in grid format). deflections are c (unit 30) written as two grids within the same file. c c (unit 31) temporary file for storing intermediate results c (unit 32) do, for harmonic continuation modes only. c c c mode 1 conversion gravity to geoid (stokes formula). c 2 gravity to deflections (vening meinesz). c 3 conversion geoid to gravity. c 4 gravity to tzz (eotvos) c 5 gravity to plumbline curvatures tyz,txz (eotvos) c 6 upward continuation to h (km) (NB! point innerzone) c c 10 downward continuation of gravity data to sea level. c 11 gravity to geoid, using harmonic continuation to a c mean height reference level, followed by stokes and c continuation of computed geoid to surface of topography. c 12 gravity to deflections using harmonic continuation. c continuation of deflections from reference level done c by using plumbline curvature parameters txz and tyz. c 0 nothing, gravity wiener filtering only c c attkm wiener filtering resolution for noisy data (kaula rule). c this option should be used only for high-pass filtering c operations (modes.ge.3), 'attkm' specifies the resolu- c tion (km) where the wiener attenuation filter is 0.5. c attkm = 0.0 means no attenuation used. c c rfi1, rla1 (degrees (or m for utm)) sw corner of wanted subgrid c if rfi1=rla1=0 then the sw-corner of the grid is used. c c inn, ine number of points of subgrid (must be even numbers) c if the wanted subgrid is bigger than the actual grid c the transform grid is padded with zeros. on output only c the non-padded part of the grid is written. c c iwndow width of cosine-tapered window zone in grid points c c c (c) rf, danish geodetic institute/national survey and cadastre c original version, june 1986 c modified and updated aug 87. c frequency domain filters changed jan 89, rf, unsw vax c last update jan 91, rf (removal of inner zone bug) c dec 93 (upw cont) c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc character gname*60,hname*60,oname*60 logical lgeo, lgeoh, lutm, lmean, ltwo, l3dec, latt dimension nnarr(2) common /gridpar/ iell,izone c c array dimensions: idim - max number of grid points c irdim - max number of points in row or column c dimension cha(2, 3080000) dimension wrk( 3200) dimension wf( 1600) idim = 3080000 iwkdim = 3200 idim2 = 2*idim c write(*,1) 1 format(/, .' ***********************************************************', .'************'/ .' * GEOFOUR - gravity field modelling by FFT - vers. jan 91 ', .'(c) RF/KMS *'/ .' ***********************************************************', .'************') c c lmean determines whether or not to remove data mean c lmean = .true. c c input file names and mode c read(*,2) gname read(*,2) hname read(*,2) oname 2 format(a60) read(*,*) mode, attkm if (mode.eq.6) read(*,*) hupw latt = (attkm.gt.0) c open(20,file=gname,status='old') if (mode.ge.10) open(21,file=hname,status='old') open(30,file=oname,status='unknown') c c unit 31 and 32 temporary scratch files c if (mode.ge.10) .open(31,status='scratch',form='unformatted') if (mode.ge.10) open(32,status='scratch' .,form='unformatted') c c constants c ltwo = (mode.eq.2.or.mode.eq.5.or.mode.eq.12) l3dec = (mode.eq.1.or.mode.eq.11) pi = 3.14159265 radeg = 180/pi degkm = 111.2 g0 = 981500.0 twodvr = 2/6371.0 c c read input data and output head c ------------------------------- c read(*,*) rfi1,rla1,inn,ine,iwndow c write(*, 9) mode, rfi1, rla1, inn, ine, iwndow 9 format(/' mode:',i3, */' wanted sw corner:',2f13.4, ', points:',2i4, *', iwndow:',i3) if (mode.ge.1.and.mode.le.6) goto (11,12,13,14,15,16),mode if (mode.ge.10) goto (17,18,19),mode-9 11 write(*,21) goto 30 12 write(*,22) goto 30 13 write(*,23) goto 30 14 write(*,24) goto 30 15 write(*,25) goto 30 16 write(*,26) hupw goto 30 17 write(*,27) goto 30 18 write(*,28) goto 30 19 write(*,29) c 21 format(' stokes formula - gravity (mgal) to geoid (meter)') 22 format(' vening-meinesz - gravity (mgal) to deflections (arcsec)') 23 format(' inverse stokes - geoid (meter) to gravity (mgal)') 24 format(' gravity (mgal) to tzz (eotvos)') 25 format(' gravity (mgal) to tyz and txz (eotvos)') 26 format(' upward continuation of gravity to ',f7.1,' km') 27 format(' downward continuation of data to sea level') 28 format(' gravity (mgal) to quasigeoid (meter) at topography') 29 format(' gravity (mgal) to deflections (") at topography') c c read data grid c -------------- c 30 write(*,36) 36 format(/' data grid:') rfic = rfi1 rlac = rla1 call rdgrid(20,rfic,rlac,inn,ine,dfi,dla,lgeo,cha,1, .ii1z,ii2z,jj1z,jj2z,idim) n = inn*ine lutm = (.not.lgeo) c c data preparation: find mean, and remove if wanted. windowing. c ------------------------------------------------------------- c r = 0.0 s = 0.0 do 31 i = 1, n r = r + cha(1,i) s = s + cha(1,i)**2 cha(2,i) = 0.0 31 continue rmean = r/n s = s/n write(*, 32) s, rmean 32 format(/' power space domain ',f10.2,', mean ',f9.2) if (lmean) write(*,34) if (.not.lmean) write(*,35) if (.not.lmean) goto 37 c do 33 i = 1, n 33 cha(1,i) = cha(1,i) - rmean 34 format(' mean value subtracted from input data prior to fft') 35 format(' mean value not removed from input data') c c windowing of data grid c 37 if (iwndow.le.0) goto 50 do 38 i = 1, iwndow 38 wf(i) = (1 - cos(3.14159265*i/(iwndow+1)))/2 wsum = 0 k = 0 sw = 0.0 do 40 i = inn, 1, -1 do 40 j = 1, ine k = k+1 w = 1.0 if (i.gt.iwndow) goto 41 w = w*wf(i) 41 if (i.le.inn-iwndow) goto 42 w = w*wf(inn-i+1) 42 if (j.gt.iwndow) goto 43 w = w*wf(j) 43 if (j.le.ine-iwndow) goto 44 w = w*wf(ine-j+1) 44 wsum = wsum + w cha(1,k) = cha(1,k)*w sw = sw + cha(1,k)**2 40 continue wsum = wsum/n sw = sw/n write(*,48) sw, wsum 48 format(' power after window ',f10.2,', wsum = ', f8.4) c c store data on unit 31 for later use c ----------------------------------- c 50 if (mode.ge.10) then do 49 i = 1, n 49 write(31) cha(1,i) endif c c fourier transformation of data c ------------------------------ c set fourier and grid constants first c internal mode parameter 'imode' controls flow of program c (continuation modes need repeat fft's - entry point at 51 or 60) c imode = mode iloop = 1 if (mode.ge.10) imode = 4 nnarr(1) = ine nnarr(2) = inn nyqn = inn/2+1 nyqe = ine/2+1 if (latt) attpi = attkm/pi if (lutm) dy = dfi*0.001 if (lutm) dx = dla*0.001 if (lgeo) dy = dfi*degkm if (lgeo) dx = dla*cos((rfic+inn/2*dfi)/radeg)*degkm c 51 if (n.le.60) call cprint(cha,idim,inn,ine) call fourt(cha,nnarr,2,-1,0,wrk,idim2,iwkdim) if (n.le.60) call cprint(cha,idim,inn,ine) c c scale transform c s = 0 do 53 i = 1,n cha(1,i) = cha(1,i)/n cha(2,i) = cha(2,i)/n s = s + cha(1,i)**2 + cha(2,i)**2 53 continue write(*, 54) s 54 format(' power freq. domain ',f10.2) c c store transform of g* for later use in mode 12 c if (mode.ne.12.or.iloop.ne.2) goto 60 rewind(31) c call setposition(31,0,0) do 56 i = 1, n 56 write(31) cha(1,i),cha(2,i) c c convolve transform with filter c ------------------------------ c c use filter including indirect effect (delta_g = -dT/dr - 2/r*T) c wiener filtering factors assume kaula rule for signal psd, c i.e. geoid data has spectral decay omega**(-4), gravity omega**(-2) c 60 du = 2*pi/dx/ine dv = 2*pi/dy/inn if (.not.latt) goto 62 omega = nyqn*dv if (mode.ne.3) write(*,61) attkm,1.0/(1+(omega*attpi)**2) if (mode.eq.3) write(*,61) attkm,1.0/(1+(omega*attpi)**4) 61 format(' wiener attenuation, resolution: ',f5.1, * ' km, att. fact. at nyquist(n): ',f6.4) c 62 do 80 iu = -nyqe+3, nyqe do 80 iv = -nyqn+3, nyqn u = (iu-1)*du v = (iv-1)*dv omega = sqrt(u**2+v**2) rr = 0.0 ri = 0.0 if (iu.eq.1.and.iv.eq.1.and.mode.ne.6) goto 79 goto (70,71,72,73,74,75,76),imode+1 c c gravity to gravity (filter only) c 70 rr = 1.0 if (latt) rr = rr/(1+(omega*attpi)**2) goto 79 c c gravity to geoid c 71 rr = 1/(omega-twodvr) if (latt) rr = rr/(1+(omega*attpi)**2) goto 79 c c deflections, pack as complex filter c 72 rr = u/(omega-twodvr) ri = v/(omega-twodvr) if (latt) rr = rr/(1+(omega*attpi)**2) if (latt) ri = ri/(1+(omega*attpi)**2) goto 79 c c geoid to gravity c 73 rr = omega-twodvr if (latt) rr = rr/(1+(omega*attpi)**4) goto 79 c c gravity to tzz c 74 rr = omega if (latt) rr = rr/(1+(omega*attpi)**2) goto 79 c c gravity to txz and tyz, complex filter c 75 rr = u ri = v if (latt) rr = rr/(1+(omega*attpi)**2) if (latt) ri = ri/(1+(omega*attpi)**2) goto 79 c c upward cont c 76 rr = exp(-omega*hupw) if (latt) rr = rr/(1+(omega*attpi)**2) c c find index and perform complex multiplication c 79 j=iu i=iv if (j.lt.1) j=j+ine if (i.lt.1) i=i+inn k = (i-1)*ine+j c c1 = cha(1,k) c2 = cha(2,k) cha(1,k) = c1*rr - c2*ri cha(2,k) = c2*rr + c1*ri 80 continue c c reverse fourier transformation c ------------------------------ c call fourt(cha,nnarr,2,1,1,wrk,idim2,iwkdim) if (n.le.60) call cprint(cha,idim,inn,ine) c c set scale constants c goto (801,81,82,83,84,84,86),imode+1 801 s = 1.0 goto 89 81 s = 1000.0/g0 goto 89 82 s = 206265/g0 goto 89 83 s = g0/1000 goto 89 84 s = 10.0 goto 89 86 s = 1.0 c 89 do 90 i = 1, n cha(1,i) = cha(1,i)*s if (ltwo) cha(2,i) = cha(2,i)*s 90 continue c c skip continuation parts in simple mode c branch to relevant action otherwise c geoid: loop 1: g->g*, loop 2: g*->geoid c defl.: loop 1: g->g*, loop 2: g*->curv corr, loop 3: g*->dfv c ------------------------------------------------------------ c if (mode.lt.10) goto 150 if (iloop.eq.1) goto 105 if (mode.eq.11) goto 120 if (mode.eq.12.and.iloop.eq.2) goto 130 if (mode.eq.12.and.iloop.eq.3) goto 135 stop c c read height data and store on unit 32 c ------------------------------------- c 105 hfic = rfi1 hlac = rla1 innh = inn ineh = ine write(*,104) 104 format(/' height grid:') call rdgrid(21,hfic,hlac,innh,ineh,dfih,dlah,lgeoh,cha,2, .ii1,ii2,jj1,jj2,idim) if (abs((dfih-dfi)/dfi).gt.0.0001.or. .abs((dlah-dla)/dla).gt.0.0001.or. .abs((rfic-hfic)/dfi).gt.0.0001.or.abs((rlac-hlac)/dla).gt. .0001.or.inn.ne.innh.or.ine.ne.ineh) .stop '*** height data grid error, must cover data points exactly' if (ii1.ne.ii1z.or.ii2.ne.ii2z.or.jj1.ne.jj1z.or.jj2.ne.jj2z) .write(*,107) 107 format(' *** warning: zero padding on gravity and height grid', .' different') c s = 0 do 109 i = 1, n s = s + cha(2,i) write(32) cha(2,i) 109 continue hmean = s/n c c continuation of gravity data to reference level c href = 0 if (mode.ne.10) href = hmean write(*, 110) href 110 format(/' harmonic continuation level: ',f8.0,' m') rewind(31) c call setposition(31,0,0) ss = 0.0 rr = 0.0 do 111 i = 1, n read(31) g if (mode.eq.10) g = g + rmean tzz = cha(1,i) h = cha(2,i) - href gstar = g + tzz*h/10000 ss = ss + (g-gstar)**2 if (abs(g-gstar).gt.rr) rr = abs(g-gstar) cha(1,i) = gstar cha(2,i) = 0 111 continue write(*,112) sqrt(ss/n), rr 112 format(' cont.corr. g to ref level: rms',f8.2,', max',f8.2, .' mgal') c c store g* and exit or jump back for new fft c if (mode.eq.10) goto 150 rewind(31) c call setposition(31, 0, 0) do 113 i = 1, n 113 write(31) cha(1,i) iloop = 2 if (mode.eq.11) imode = 1 if (mode.eq.12) imode = 5 goto 51 c c geoid: continue computed geoid to surface of topography c 120 continue rewind(31) rewind(32) c call setposition(31,0,0) c call setposition(32,0,0) rr = 0.0 ss = 0.0 do 121 i = 1, n read(31) g read(32) h r = -(g+rmean)*(h-href)/g0 cha(1,i) = cha(1,i) + r ss = ss + r**2 if (rr.lt.abs(r)) rr = abs(r) 121 continue write(*,122) sqrt(ss/n),rr 122 format(' cont.coor. geoid to topo: rms',f8.3,', max',f8.3 .,' m') goto 150 c c deflections - loop 2: store curvature corrections on unit 32 c 130 s = 206265.0/g0/10000 rewind(31) rewind(32) c call setposition(31, 0, 0) c call setposition(32, 0, 0) do 131 i = 1, n read(32) h cha(1,i) = -cha(1,i)*(h-href)*s cha(2,i) = -cha(2,i)*(h-href)*s 131 continue rewind(32) c call setposition(32, 0, 0) do 132 i = 1, n 132 write(32) cha(1,i),cha(2,i) c iloop = 3 imode = 2 do 133 i = 1, n 133 read(31) cha(1,i),cha(2,i) goto 60 c c deflections - loop 3: add curvature corrections to dfv result c 135 continue rewind(32) c call setposition(32,0,0) ss = 0.0 rr = 0.0 do 137 i = 1, n read(32) r,s cha(1,i) = cha(1,i) + r cha(2,i) = cha(2,i) + s ss = ss + r**2 + s**2 if (abs(r).gt.rr) rr = abs(r) if (abs(s).gt.rr) rr = abs(s) 137 continue write(*,138) sqrt(ss/(2*n)), rr 138 format(' cont.corr. dfv to topo: rms',f8.2,', max',f8.2, .' arcsec') goto 150 c c find min and max of computed values c ----------------------------------- c 150 nk = 1 if (ltwo) nk = 2 do 155 k = 1, nk c1 = 999999.9 c2 = -999999.9 do 151 i = 1, n if (cha(k,i).lt.c1) c1 = cha(k,i) if (cha(k,i).gt.c2) c2 = cha(k,i) 151 continue write(*,152) c1, c2 152 format(/' min and max computed values: ',2f11.2) 155 continue c c write out results on ofile - without possible zero padding c ---------------------------------------------------------- c do 920 k = 1, nk if (lgeo) write(30,901) rfic+ii1z*dfi,rfic+(inn-1-ii2z)*dfi, *rlac+jj1z*dla,rlac+(ine-1-jj2z)*dla,dfi,dla if (lutm) write(30,902) rfic+ii1z*dfi,rfic+(inn-1-ii2z)*dfi, *rlac+jj1z*dla,rlac+(ine-1-jj2z)*dla,dfi,dla,iell,izone 901 format(/' ',4f12.6,2f12.7) 902 format(/' ',6f11.0,/,' ',2i5) do 910 i = 1+ii2z, inn-ii1z if (.not.l3dec) . write(30,912) (cha(k,(i-1)*ine+j),j=1+jj1z,ine-jj2z) if (l3dec) . write(30,913) (cha(k,(i-1)*ine+j),j=1,ine-jj2z) 910 continue 912 format(/,40(' ',8f9.2,/)) 913 format(/,40(' ',8f9.3,/)) 920 continue c close(31,status='delete') if (mode.gt.10) close(32,status='delete') end c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c p r i n t c c prints the contents of complex array cha c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine cprint(cha, ihadim, inn, ine) implicit double precision(a-h,o-z) dimension cha(2,ihadim) do 19 k=1,2 j1=1 j2=ine write(*,11) 11 format(' ') do 18 i=inn,1,-1 write(*,12) (cha(k,j),j=j1,j2) 12 format(8f11.4) j1 = j1+ine j2 = j2+ine 18 continue 19 continue return end c subroutine fourt(datt,nn,ndim,isign,iform,work, .idim1,idim2) implicit double precision(a-h,o-z) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c f o u r t c c version=740301 c program description norsar n-pd9 dated 1 july 1970 c author n m brenner›A c further description three fortran programs etc. c issued by lincoln laboratory, mit, july 1967 c two corrections by hjortenberg 1974 c the fast fourier transform in usasi basic fortran c c modified to rc fortran rf june 84 c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc dimension datt(idim1),nn(ndim),ifact(32),work(idim2) c np0=0 nprev=0 c twopi=6.283185307d0 rthlf=.7071067812d0 if(ndim-1)920,1,1 1 ntot=2 do 2 idim=1,ndim if(nn(idim))920,920,2 2 ntot=ntot*nn(idim) c c mainloop for each dimension c np1=2 do 910 idim=1,ndim n=nn(idim) np2=np1*n if(n-1)920,900,5 c c is n a power of two and if not, what are its factors c 5 m=n ntwo=np1 if=1 idiv=2 10 iquot=m/idiv irem=m-idiv*iquot if(iquot-idiv)50,11,11 11 if(irem)20,12,20 12 ntwo=ntwo+ntwo ifact(if)=idiv if=if+1 m=iquot go to 10 20 idiv=3 inon2=if 30 iquot=m/idiv irem=m-idiv*iquot if(iquot-idiv)60,31,31 31 if(irem)40,32,40 32 ifact(if)=idiv if=if+1 m=iquot go to 30 40 idiv=idiv+2 go to 30 50 inon2=if if(irem)60,51,60 51 ntwo=ntwo+ntwo go to 70 60 ifact(if)=m 70 non2p=np2/ntwo c c separate four cases--- c 1. complex transform c 2. real transform for the 2nd, 3nd, etc. dimension. method-- c transform half the datt, supplying the other half by con- c jugate symmetry. c 3. real transform for the 1st dimension,n odd. method-- c set the imaginary parts to zero c 4. real transform for the 1st dimension,n even.method-- c transform a complex array of lenght n/2 whose real parts c are the even numberd real values and whose imaginary parts c are the odd numberedreal values. separate and supply c the second half by conjugate summetry. c icase=1 ifmin=1 i1rng=np1 if(idim-4)74,100,100 74 if(iform)71,71,100 71 icase=2 i1rng=np0*(1+nprev/2) if(idim-1)72,72,100 72 icase=3 i1rng=np1 if(ntwo-np1)100,100,73 73 icase=4 ifmin=2 ntwo=ntwo/2 n=n/2 np2=np2/2 ntot=ntot/2 i=1 do 80 j=1,ntot datt(j)=datt(i) 80 i=i+2 c c shuffle datt by bit reversal, since n=2**k. as the shuffling c can be done by simple interchange, no working array is needed c 100 if(non2p-1)101,101,200 101 np2hf=np2/2 j=1 do 150 i2=1,np2,np1 if(j-i2)121,130,130 121 i1max=i2+np1-2 do 125 i1=i2,i1max,2 do 125 i3=i1,ntot,np2 j3=j+i3-i2 tempr=datt(i3) tempi=datt(i3+1) datt(i3)=datt(j3) datt(i3+1)=datt(j3+1) datt(j3)=tempr 125 datt(j3+1)=tempi 130 m=np2hf 140 if(j-m)150,150,141 141 j=j-m m=m/2 if(m-np1)150,140,140 150 j=j+m go to 300 c c shuffle datt by digit reversal for general n c 200 nwork=2*n do 270 i1=1,np1,2 do 270 i3=i1,ntot,np2 j=i3 do 260 i=1,nwork,2 if(icase-3)210,220,210 210 work(i)=datt(j) work(i+1)=datt(j+1) go to 240 220 work(i)=datt(j) work(i+1)=0. 240 ifp2=np2 if=ifmin 250 ifp1=ifp2/ifact(if) j=j+ifp1 if(j-i3-ifp2)260,255,255 255 j=j-ifp2 ifp2=ifp1 if=if+1 if(ifp2-np1)260,260,250 260 continue i2max=i3+np2-np1 i=1 do 270 i2=i3,i2max,np1 datt(i2)=work(i) datt(i2+1)=work(i+1) 270 i=i+2 c c main loop for factors of two c w=exp(isign*2*pi*sqrt(-1)*m/(4*mmax)). check for w=isign*sqrt(-1) c and repeat for w=w*(1+isign*sqrt(-1))/sqrt(2) c 300 if(ntwo-np1)600,600,305 305 np1tw=np1+np1 ipar=ntwo/np1 310 if(ipar-2)350,330,320 320 ipar=ipar/4 go to 310 330 do 340 i1=1,i1rng,2 do 340 k1=i1,ntot,np1tw k2=k1+np1 tempr=datt(k2) tempi=datt(k2+1) datt(k2)=datt(k1)-tempr datt(k2+1)=datt(k1+1)-tempi datt(k1)=datt(k1)+tempr 340 datt(k1+1)=datt(k1+1)+tempi 350 mmax=np1 360 if(mmax-ntwo/2)370,600,600 370 lmax=max0(np1tw,mmax/2) do 570 l=np1,lmax,np1tw m=l if(mmax-np1)420,420,380 380 theta=-twopi*dble(l)/dble(4*mmax) if(isign)400,390,390 390 theta=-theta 400 wr=cos(theta) wi=sin(theta) 410 w2r=wr*wr-wi*wi w2i=2.*wr*wi w3r=w2r*wr-w2i*wi w3i=w2r*wi+w2i*wr 420 do 530 i1=1,i1rng,2 kmin=i1+ipar*m if(mmax-np1)430,430,440 430 kmin=i1 440 kdif=ipar*mmax 450 kstep=4*kdif if(kstep-ntwo)460,460,530 460 do 520 k1=kmin,ntot,kstep k2=k1+kdif k3=k2+kdif k4=k3+kdif if(mmax-np1)470,470,480 470 u1r=datt(k1)+datt(k2) u1i=datt(k1+1)+datt(k2+1) u2r=datt(k3)+datt(k4) u2i=datt(k3+1)+datt(k4+1) u3r=datt(k1)-datt(k2) u3i=datt(k1+1)-datt(k2+1) if(isign)471,472,472 471 u4r=datt(k3+1)-datt(k4+1) u4i=datt(k4)-datt(k3) go to 510 472 u4r=datt(k4+1)-datt(k3+1) u4i=datt(k3)-datt(k4) go to 510 480 t2r=w2r*datt(k2)-w2i*datt(k2+1) t2i=w2r*datt(k2+1)+w2i*datt(k2) t3r=wr*datt(k3)-wi*datt(k3+1) t3i=wr*datt(k3+1)+wi*datt(k3) t4r=w3r*datt(k4)-w3i*datt(k4+1) t4i=w3r*datt(k4+1)+w3i*datt(k4) u1r=datt(k1)+t2r u1i=datt(k1+1)+t2i u2r=t3r+t4r u2i=t3i+t4i u3r=datt(k1)-t2r u3i=datt(k1+1)-t2i if(isign)490,500,500 490 u4r=t3i-t4i u4i=t4r-t3r go to 510 500 u4r=t4i-t3i u4i=t3r-t4r 510 datt(k1)=u1r+u2r datt(k1+1)=u1i+u2i datt(k2)=u3r+u4r datt(k2+1)=u3i+u4i datt(k3)=u1r-u2r datt(k3+1)=u1i-u2i datt(k4)=u3r-u4r 520 datt(k4+1)=u3i-u4i kdif=kstep kmin=4*(kmin-i1)+i1 go to 450 530 continue m=m+lmax if(m-mmax)540,540,570 540 if(isign)550,560,560 550 tempr=wr wr=(wr+wi)*rthlf wi=(wi-tempr)*rthlf go to 410 560 tempr=wr wr=(wr-wi)*rthlf wi=(tempr+wi)*rthlf go to 410 570 continue ipar=3-ipar mmax=mmax+mmax go to 360 c c main loop for factoers not equal to two c w=exp(isign*2*pi*sqrt(-1)*(j1+j2-i3-1)/ifp2) c 600 if(non2p-1)700,700,601 601 ifp1=ntwo if=inon2 610 ifp2=ifact(if)*ifp1 theta=-twopi/dble(ifact(if)) if(isign)612,611,611 611 theta=-theta 612 wstpr=cos(theta) wstpi=sin(theta) do 650 j1=1,ifp1,np1 thetm=-twopi*dble(j1-1)/dble(ifp2) if(isign)614,613,613 613 thetm=-thetm 614 wminr=cos(thetm) wmini=sin(thetm) i1max=j1+i1rng-2 do 650 i1=j1,i1max,2 do 650 i3=i1,ntot,np2 i=1 wr=wminr wi=wmini j2max=i3+ifp2-ifp1 do 640 j2=i3,j2max,ifp1 twowr=wr+wr j3max=j2+np2-ifp2 do 630 j3=j2,j3max,ifp2 jmin=j3-j2+i3 j=jmin+ifp2-ifp1 sr=datt(j) si=datt(j+1) oldsr=0. oldsi=0. j=j-ifp1 620 stmpr=sr stmpi=si sr=twowr*sr-oldsr+datt(j) si=twowr*si-oldsi+datt(j+1) oldsr=stmpr oldsi=stmpi j=j-ifp1 if(j-jmin)621,621,620 621 work(i)=wr*sr-wi*si-oldsr+datt(j) work(i+1)=wi*sr+wr*si-oldsi+datt(j+1) 630 i=i+2 wtemp=wr*wstpi wr=wr*wstpr-wi*wstpi 640 wi=wi*wstpr+wtemp i=1 do 650 j2=i3,j2max,ifp1 j3max=j2+np2-ifp2 do 650 j3=j2,j3max,ifp2 datt(j3)=work(i) datt(j3+1)=work(i+1) 650 i=i+2 if=if+1 ifp1=ifp2 if(ifp1-np2)610,700,700 c c complete areal transform in the 1st dimension, n even, by con- c jugate symmetries c 700 go to (900,800,900,701),icase 701 nhalf=n n=n+n theta=-twopi/dble(n) if(isign)703,702,702 702 theta=-theta 703 wstpr=cos(theta) wstpi=sin(theta) wr=wstpr wi=wstpi imin=3 jmin=2*nhalf-1 go to 725 710 j=jmin do 720 i=imin,ntot,np2 sumr=(datt(i)+datt(j))/2. sumi=(datt(i+1)+datt(j+1))/2. difr=(datt(i)-datt(j))/2. difi=(datt(i+1)-datt(j+1))/2. tempr=wr*sumi+wi*difr tempi=wi*sumi-wr*difr datt(i)=sumr+tempr datt(i+1)=difi+tempi datt(j)=sumr-tempr datt(j+1)=-difi+tempi 720 j=j+np2 imin=imin+2 jmin=jmin-2 wtemp=wr*wstpi wr=wr*wstpr-wi*wstpi wi=wi*wstpr+wtemp 725 if(imin-jmin)710,730,740 730 if(isign)731,740,740 731 do 735 i=imin,ntot,np2 735 datt(i+1)=-datt(i+1) 740 np2=np2+np2 ntot=ntot+ntot j=ntot+1 imax=ntot/2+1 745 imin=imax-2*nhalf i=imin go to 755 750 datt(j)=datt(i) datt(j+1)=-datt(i+1) 755 i=i+2 j=j-2 if(i-imax)750,760,760 760 datt(j)=datt(imin)-datt(imin+1) datt(j+1)=0. if(i-j)770,780,780 765 datt(j)=datt(i) datt(j+1)=datt(i+1) 770 i=i-2 j=j-2 if(i-imin)775,775,765 775 datt(j)=datt(imin)+datt(imin+1) datt(j+1)=0. imax=imin go to 745 780 datt(1)=datt(1)+datt(2) datt(2)=0. go to 900 c c complete a real transform for the 2nd, 3rd, etc. dimension by c conjugate symmetries. c 800 if(i1rng-np1)805,900,900 805 do 860 i3=1,ntot,np2 i2max=i3+np2-np1 do 860 i2=i3,i2max,np1 imax=i2+np1-2 imin=i2+i1rng jmax=2*i3+np1-imin if(i2-i3)820,820,810 810 jmax=jmax+np2 820 if(idim-2)850,850,830 830 j=jmax+np0 do 840 i=imin,imax,2 datt(i)=datt(j) datt(i+1)=-datt(j+1) 840 j=j-2 850 j=jmax do 860 i=imin,imax,np0 datt(i)=datt(j) datt(i+1)=-datt(j+1) 860 j=j-np0 c c end of loop on each dimension c 900 np0=np1 np1=np2 910 nprev=n 920 return end c subroutine rdgrid(iunit, rfic, rlac, inn, ine, dfi, dla, .lgeo, cha, ik, ii1z, ii2z, jj1z, jj2z, idim) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c r d g r i d c c subroutine for reading a digital grid file on c standard format, i.e. stored rowwise from nw to se, with label. c c as se-corner coordinate 'rfic, rlac' (degrees) will be used c (unless they are zero, then the grid corner is used). c a grid containing 'inn' x 'ine' points will be put in array c 'cha' of declared dimension 'idim'. c if inn=0 the complete grid will be read. c if the wanted grid is too large a zero padding will be done. c c last updated jun 90, rf c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc implicit double precision(a-h,o-z) common /gridpar/ iell,izone dimension cha(2, idim) dimension hlab(6), hrow(3200) logical lgeo, lutm irdim = 3200 c c initialize statistics c nr = 0 rsum = 0.0 rsum2 = 0.0 rmin = 99999 rmax = -99999 c read(iunit,*) (hlab(j),j=1,6) lutm = (abs(hlab(1)).ge.100.or.abs(hlab(2)).ge.100) lgeo = (.not.lutm) if (lgeo) goto 111 read(iunit,*) iell,izone write(*,110) iell,izone if (iell.lt.1.or.iell.gt.3.or.izone.lt.1.or.izone.gt.60) *stop 'illegal ellipsoid or utm zone' 110 format(' - input grid in utm, ell ',i1,' zone ',i2,' -') 111 dfi = hlab(5) dla = hlab(6) nn = (hlab(2)-hlab(1))/dfi+1.5 ne = (hlab(4)-hlab(3))/dla+1.5 if (nn.gt.irdim.or.ne.gt.irdim) stop 'too long rows/columns' c c find corner indices for wanted subgrid c if (inn.eq.0) then inn = nn ine = ne endif if (rfic.eq.0.and.rlac.eq.0) then rfic = hlab(1) rlac = hlab(3) endif ifi1 = (rfic-hlab(1))/dfi+1.5 ila1 = (rlac-hlab(3))/dla+1.5 ifi2 = ifi1+inn-1 ila2 = ila1+ine-1 rfic = (ifi1-1)*dfi + hlab(1) rlac = (ila1-1)*dla + hlab(3) n = inn*ine c c check boundaries for padding c ii1z = 0 ii2z = 0 jj1z = 0 jj2z = 0 if (ifi1.lt.1) ii1z = 1-ifi1 if (ifi2.gt.nn) ii2z = ifi2-nn if (ila1.lt.1) jj1z = 1-ila1 if (ila2.gt.ne) jj2z = ila2-ne c if (n.gt.idim) then write(*, 122) n,idim 122 format(' *** array dim too small - wanted, declared ',2i8) stop ' *** sorry ***' endif c c read data grid values c data in cha array stored with first element at nw corner c ilast = ifi1 if (ilast.lt.1) ilast = 1 c do 130 i = nn,ilast,-1 c read(iunit,*,end=131) (hrow(j),j=1,ne) c if (i.lt.ifi1.or.i.gt.ifi2) goto 130 jj0 = (ifi2-i)*ine - ila1+1 do 129 j = 1,ne r = hrow(j) if (j.lt.ila1.or.j.gt.ila2) goto 129 cha(ik,j+jj0) = r nr = nr + 1 if (r.gt.rmax) rmax = r if (r.lt.rmin) rmin = r rsum = rsum + r rsum2 = rsum2 + r**2 129 continue 130 continue goto 133 131 write(*,132) i 132 format(' *** too few data in grid file, lastrow = ',i7) stop ' *** check grid label and data ***' c c zero padding c 133 if (ii1z+ii2z+jj1z+jj2z.gt.0) then do 138 i = inn,1,-1 jj0 = (inn-i)*ine if (i.gt.ifi2.or.i.lt.ifi1) then do 134 j = 1, ine 134 cha(ik,j+jj0) = 0 else do 135 j = 1, jj1z 135 cha(ik,j+jj0) = 0 do 136 j = ine-jj2z+1, ine 136 cha(ik,j+jj0) = 0 endif 138 continue endif c c write information and statistics c if (nr.eq.0) stop '*** no points read from grid, wrong area' rfi = hlab(1) + (ifi1-1)*dfi rla = hlab(3) + (ila1-1)*dla r = rsum/nr s = 0.0 if (n.gt.1) .s = sqrt((rsum2 - rsum**2/nr)/(nr-1)) if (lgeo) write(*,141) (hlab(j),j=1,6),nn,ne if (lutm) write(*,142) (hlab(j),j=1,6),nn,ne 141 format(' gridlab:',4f10.4,2f9.4,i5,i4) 142 format(' gridlab:',4f10.0,2f8.0,i5,i4) if (lgeo) write(*, 143) rfi, rla, inn, ine, inn*ine 143 format(' selected: sw corner ',2f10.4, ', points ', 2i6,i8) if (lutm) write(*, 144) rfi, rla, inn, ine, ii 144 format(' selected: sw corner ',2f10.0, ', points ', 2i6,i8) write(*, 145) nr, r, s, rmin, rmax 145 format(' statistics of data selected from input grid:' ./' pts mean std.dev. min max:',i8,4f9.2) if (ii1z+ii2z+jj1z+jj2z.gt.0) write(*,146) ii1z,ii2z,jj1z,jj2z 146 format(' zero padding done on grid, no of rows/cols S/N/E/W:',4i4) return end