c--------( RPDEQS : calc. reduction-to-pole by source )-------- parameter (LHW=100, LXO=1201, LYO=1201) parameter (LXM=LXO+LHW+LHW, LYM=LYO+LHW+LHW) parameter (DPI=3.14159*2.) c dimension s(LXM,LYM), t(LXM,LYM), f(LXO,LYO), h(LXO,LYO) character*80 wdr, flog, fnam1, fnam2, fnam3 character buf*80, out*72, area*8, area1*8, area2*8, sdtm*20 c call premsg('RPDEQS : calc. reduc-to-pole by source') call opnpin() ldr = lwkdir(50,wdr) + 1 flog = wdr fnam1 = wdr fnam2 = wdr fnam3 = wdr call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80)) call gparma('Enter Reduction Alt. filename ==> ', * 81-ldr, fnam1(ldr:80)) open(1,file=fnam1,status='old') 1 read(1,'(a)') buf if(buf(1:1).eq.'#') goto 1 read(buf,'(a8,i4,4x,4i8)') area, ncc, iorg,korg, ispa,ispb nc = ncc if(ncc.ge.200) nc = ncc - 200 if(nc.ge.1.and.nc.le.62) then iorg = 0 ispa = 0 ispb = 0 korg = (nc-30)*360 - 180 if(nc.gt.60) korg = 0 if(nc.eq.61) iorg = +5400 if(nc.eq.62) iorg = -5400 endif read(1,'(a)') buf read(buf,*) ixl, iyl, mszx, mszy, mx, my, hnul if(mx.gt.LXO.or.my.gt.LYO) call abendm('too large array size') if(mx.le.0.or.my.le.0) call abendm('illegal array size') read(1,*) ((h(i,j),i=1,mx),j=1,my) close(1) xl = float(ixl)/1000. yl = float(iyl)/1000. ixh = ixl + mszx*(mx-1) iyh = iyl + mszy*(my-1) c call gparma('Enter ADEQS model input filename ==> ', * 81-ldr, fnam2(ldr:80)) open(2,file=fnam2,status='old') 2 read(2,'(a)') buf if(buf(1:1).eq.'#') goto 2 read(buf,'(a8,i4,4x,4i8)') area1, ncc1, iorg1,korg1, ispa1,ispb1 nc1 = ncc1 if(mcc1.ge.200) nc1 = ncc1 - 200 if(nc1.ge.1.and.nc1.le.62) then iorg1 = 0 ispa1 = 0 ispb1 = 0 korg1 = (nc1-30)*360 - 180 if(nc1.gt.60) korg1 = 0 if(nc1.eq.61) iorg1 = +5400 if(nc1.eq.62) iorg1 = -5400 endif read(2,*) ixlm, iylm, mszx1, mszy1, mxm, mym, anul, dmy, klp if(ncc1.ne.ncc .or. * iorg1.ne.iorg.or.korg1.ne.korg .or. * ispa1.ne.ispa.or.ispb1.ne.ispb .or. * mszx1.ne.mszx.or.mszy1.ne.mszy) * call abendm('Header inconsistent') if(mxm.gt.LXM.or.mym.gt.LYM) call abendm('too large array size') if(mxm.lt.10.or.mym.lt.10) call abendm('illegal array size') read(2,*) ((s(i,j),i=1,mxm),j=1,mym) 3 read(2,'(a)') buf if(buf(1:1).eq.'#') goto 3 read(buf,'(a8,i4,4x,4i8)') area2, ncc1, iorg1,korg1, ispa1,ispb1 nc1 = ncc1 if(ncc1.ge.200) nc1 = ncc1 - 200 if(nc1.ge.1.and.nc1.le.62) then iorg1 = 0 ispa1 = 0 ispb1 = 0 korg1 = (nc1-30)*360 - 180 if(nc1.gt.60) korg1 = 0 if(nc1.eq.61) iorg1 = +5400 if(nc1.eq.62) iorg1 = -5400 endif read(2,'(a)') buf read(buf,*) ixl1, iyl1, mszx1, mszy1, mxm1, mym1, tnul if(area2.ne.area1.or.ncc1.ne.ncc .or. * iorg1.ne.iorg.or.korg1.ne.korg .or. * ispa1.ne.ispa.or.ispb1.ne.ispb) * call abendm('Header inconsistent') if(ixl1.ne.ixlm.or.iyl1.ne.iylm .or. * mszx1.ne.mszx.or.mszy1.ne.mszy .or. * mxm1.ne.mxm.or.mym1.ne.mym) * call abendm('Header inconsistent') read(2,*) ((t(i,j),i=1,mxm),j=1,mym) read(2,'(2x,4f8.1)') dip, dec, the, phi close(2) xlm = float(ixlm)/1000. ylm = float(iylm)/1000. ixhm = ixlm + mszx*(mxm-1) iyhm = iylm + mszy*(mym-1) if(ixl.le.ixlm.or.ixh.ge.ixhm .or. * iyl.le.iylm.or.iyh.ge.iyhm) call abendm('ranges conflict') smx = float(mszx) smy = float(mszy) amxy = smx * smy vol = amxy * sqrt(amxy) c write(out,'(a,a8,a,i4,4x,4i8)') * '--- Area : ', area, ' Coord.', ncc, iorg,korg, ispa,ispb call premsg(out) write(out,'(a,2a10,2x,2a8,2a6,a7)') * '--- (mesh) ', 'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my', * 'alt' call premsg(out) write(out,'(a,2f10.2,2i8,2i6)') * '--- reduc. alt. ', xl, yl, mszx,mszy, mx, my call premsg(out) write(out,'(a,2f10.2,2i8,2i6,a8)') * '--- equiv.source', xlm,ylm, mszx,mszy, mxm,mym, ' var.' call premsg(out) call premsg(' Truncation of Source effect ?') call gparmf(' Half width (in km) : ', hw) hwm = hw * 1000. nhx = nint(hwm/smx) nhy = nint(hwm/smy) c call gparma('Enter Reduc-to-Pole out filename ==> ', * 81-ldr, fnam3(ldr:80)) call clspin() call strdtm(sdtm) open(3,file=fnam3,status='new') write(6,'(a)') '----------------------------------------' write(6,'(a)') 'RPDEQS : calc. reduc-to-pole by source' write(6,'(a,a)') * ' Reduction alt. infile : ', fnam1(1:lrtrim(fnam1)) write(6,'(a,a)') * ' ADEQS model infile : ', fnam2(1:lrtrim(fnam2)) write(6,'(a,a)') * ' Reduc-to-Pole outfile : ', fnam3(1:lrtrim(fnam3)) write(6,1000) area, ncc, iorg,korg, ispa,ispb, * xl, yl, mszx, mszy, mx, my, * xlm, ylm, mszx, mszy, mxm, mym, ' var.' write(6,2000) hw, nhx, nhy c call dpcini('... calc. Reduc-Pole anom.: ') do 200 j=1,my do 200 i=1,mx if(h(i,j).eq.hnul) then f(i,j) = 9999.9 else xp = float(ixl - ixlm + mszx*(i-1)) yp = float(iyl - iylm + mszy*(j-1)) ib = nint(xp/smx) + 1 jb = nint(yp/smy) + 1 ss = 0. do 100 kj=-nhy,+nhy js = jb + kj y = float(kj)*smy do 100 ki=-nhx,+nhx is = ib + ki x = float(ki)*smx if(is.le.0.or.is.gt.mxm) goto 100 if(js.le.0.or.js.gt.mym) goto 100 if(t(is,js).eq.tnul) goto 100 z = h(i,j) - t(is,js) rr = x*x + y*y + z*z r = sqrt(rr) v = z*z*2. - x*x - y*y ss = ss + s(is,js)*v*vol/(rr*rr*r) 100 continue f(i,j) = ss endif call dpcent((j-1)*mx+i, mx*my) 200 continue c write(3,'(a8,i4,4x,4i8)') area, ncc, iorg,korg, ispa,ispb write(3,'(2i12,4i6,1x,f7.1,1x,f7.0)') * ixl, iyl, mszx, mszy, mx, my, 9999.9, 0. do 300 j=1,my write(3,'((f7.1,9(1x,f7.1)))') (f(i,j),i=1,mx) 300 continue write(3,'(a8,i4,4x,4i8)') area, ncc, iorg,korg, ispa,ispb write(3,'(2i12,4i6,1x,f7.1,1x,f7.0)') * ixl, iyl, mszx, mszy, mx, my, hnul, -1. do 400 j=1,my write(3,'((f7.1,9(1x,f7.1)))') (h(i,j),i=1,mx) 400 continue close(3) c if(flog(ldr:ldr).ne.' ') then open(99,file=flog,access='append') write(99,'(//a)') * 'RPDEQS : reduc-to-pole anomaly by source' write(99,'(a,a)') * ' Reduction alt. infile : ', fnam1(1:lrtrim(fnam1)) write(99,'(a,a)') * ' ADEQS model infile : ', fnam2(1:lrtrim(fnam2)) write(99,'(a,a)') * ' Reduc-to-Pole outfile : ', fnam3(1:lrtrim(fnam3)) write(99,'(a,a,a)') '--------------- ', sdtm, '---------------' write(99,1000) area, ncc, iorg,korg, ispa,ispb, * xl, yl, mszx, mszy, mx, my, * xlm, ylm, mszx, mszy, mxm, mym, ' var.' write(99,2000) hw, nhx, nhy call strdtm(sdtm) write(99,'(a,a,a)') '=============== ', sdtm, '===============' close(99) endif stop c 1000 format('Areaname : ', a8, 5x, 'Coord.', i4, 4x, 4i8/ * 3x, ' (mesh) ', ' xs ', ' ys ', * ' msz-x', ' msz-y', ' mx', ' my', ' alt '/ * 3x, 'reduc. surf.', 2f10.2, 2i8, 2i6/ * 3x, 'equiv.source', 2f10.2, 2i8, 2i6, a8) 2000 format(' Half width of calc.window :', f6.2, * ' km (', i4, ' *', i4, ' source-mesh )') end