c--------( CMAGF : calc. anomaly of fine model )-------- 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), e(6) character*80 wdr, flog, fnam1, fnam2, fnam3, fnam4 character buf*80, out*72, area*8, area1*8, area2*8 character strs*20, sdtm*20 c call premsg('CMAGF : calc. anomaly of fine model') call opnpin() ldr = lwkdir(50,wdr) + 1 flog = wdr fnam1 = wdr fnam2 = wdr fnam3 = wdr fnam4 = wdr call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80)) c call gparma('Enter Calc.Alt. data 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,*) ixs, iys, 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') xs = float(ixs)/1000. ys = float(iys)/1000. ixt = ixs + mszx*(mx-1) iyt = iys + mszy*(my-1) read(1,*) ((h(i,j),i=1,mx),j=1,my) close(1) smx = float(mszx) smy = float(mszy) c call gparma('Enter AMAG 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(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,*) ixl, iyl, mszxm, mszym, 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) * 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') xl = float(ixl)/1000. yl = float(iyl)/1000. ixh = ixl + mszxm*(mxm-1) iyh = iyl + mszym*(mym-1) if(ixl.gt.ixs.or.iyl.gt.iys .or. * ixh.lt.ixt.or.iyh.lt.iyt) call abendm('ranges conflict') read(2,*) ((s(i,j),i=1,mxm),j=1,mym) close(2) smxm = float(mszxm) smym = float(mszym) c call gparma('Enter Source Alt. data filename ==> ', * 81-ldr, fnam3(ldr:80)) open(10,file=fnam3,status='old') 3 read(10,'(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(10,'(a)') buf read(buf,*) ixst, iyst, mszxt, mszyt, mxt, myt, tnul if(ncc1.ne.ncc .or. * iorg1.ne.iorg.or.korg1.ne.korg .or. * ispa1.ne.ispa.or.ispb1.ne.ispb) * call abendm('Header inconsistent') if(mxt.gt.LXM.or.myt.gt.LYM) call abendm('too large array size') if(mxt.lt.10.or.myt.lt.10) call abendm('illegal array size') nmsx = mszxm / mszxt if((nmsx*mszxt).ne.mszxm) call abendm('mesh size illegal') nmsy = mszym / mszyt if((nmsy*mszyt).ne.mszym) call abendm('mesh size illegal') hnmsx = float(nmsx-1)/2. hnmsy = float(nmsy-1)/2. xst = float(ixst)/1000. yst = float(iyst)/1000. ixtt = ixst + mszxt*(mxt-1) iytt = iyst + mszyt*(myt-1) if(ixst.gt.ixl.or.iyst.gt.iyl .or. * ixtt.lt.ixh.or.iytt.lt.iyh) call abendm('ranges conflict') read(10,*) ((t(i,j),i=1,mxt),j=1,myt) close(10) smxt = float(mszxt) smyt = float(mszyt) 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)') * '--- (mesh) ', 'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my' call premsg(out) write(out,'(a,2f10.2,2i8,2i6)') * '--- Calc. Alt. ', xs, ys, mszx,mszy, mx,my call premsg(out) write(out,'(a,2f10.2,2i8,2i6,a8)') * '--- AMAG model ', xl, yl, mszxm,mszym, mxm,mym write(out,'(a,2f10.2,2i8,2i6)') * '--- Source Alt. ', xst, yst, mszxt,mszyt, mxt,myt 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 premsg('') call premsg('Specify Source-Bottom configuration') call prompt(' Constant Thickness (1) or Flat Bottom (0) ') call gparmi(': ', ks) if(ks.lt.0.or.ks.gt.1) call abendm('invalid value') strs = ' below sea level ' if(ks.eq.1) strs = ' below terr.surface ' call gparmf(' Bottom-depth (in m)' // strs // ': ', botm) call premsg('Specify Direction Parameters (in degrees)') call gparmf2(' Inclin., Declin. of Ambient field : ', dip, dec) if(dip.lt.-90..or.dip.gt.90. .or. * dec.lt.-180..or.dec.ge.360.) call abendm('invalid value') call gparmf2(' Inclin., Declin. of Magnetization : ', the, phi) if(the.lt.-90..or.the.gt.90. .or. * phi.lt.-180..or.phi.ge.360.) call abendm('invalid value') c call gparma('Enter Calc.Anom. output filename ==> ', * 81-ldr, fnam4(ldr:80)) call clspin() call strdtm(sdtm) open(3,file=fnam4,status='new') write(6,'(a)') '----------------------------------------' write(6,'(a)') 'CMAGF : calc. anomaly of fine model' write(6,'(a,a)') * ' Calc.Alt. data infile : ', fnam1(1:lrtrim(fnam1)) write(6,'(a,a)') * ' AMAG model infile : ', fnam2(1:lrtrim(fnam2)) write(6,'(a,a)') * ' Source Alt. infile : ', fnam3(1:lrtrim(fnam3)) write(6,'(a,a)') * ' Calc.Anom. outfile : ', fnam4(1:lrtrim(fnam4)) write(6,1000) area, ncc, iorg,korg, ispa,ispb, * xs, ys, mszx, mszy, mx, my, * xl, yl, mszxm,mszym, mxm,mym, * xst,yst, mszxt,mszyt, mxt,myt write(6,2000) hw, nhx, nhy c do 100 j=1,my do 100 i=1,mx f(i,j) = 0. if(h(i,j).eq.hnul) f(i,j) = 9999.9 100 continue call dpcini('... calc. model anomaly : ') call magafd(the, phi, dip, dec) am = smxt * smyt do 200 j=1,mym iyc = iyl + (j-1)*mszym ja = nint(float(iyc-iyst)/smyt) - hnmsy kjc = nint(float(iyc-iys)/smy) do 200 i=1,mxm ixc = ixl + (i-1)*mszxm ia = nint(float(ixc-ixst)/smxt) - hnmsx kic = nint(float(ixc-ixs)/smx) do 300 js=1,nmsy jt = ja + js if(jt.le.0.or.jt.gt.myt) goto 300 iyt = iyst + (jt-1)*mszyt do 400 is=1,nmsx it = ia + is if(it.le.0.or.it.gt.mxt) goto 400 ixt = ixst + (it-1)*mszxt ht = t(ia+1,ja+1) if(ht.eq.tnul) call abendm('source alt. undefined') do 500 kj=-nhy,nhy jr = kjc + kj if(jr.lt.0.or.jr.ge.my) goto 500 y = float(iyt - iys - jr*mszy) do 600 ki=-nhx,nhx ir = kic + ki x = float(ixt - ixs - ir*mszx) if(ir.lt.0.or.ir.ge.mx) goto 600 hc = h(ir+1,jr+1) if(hc.eq.hnul) goto 600 za = hc - ht sz = botm if(ks.eq.0) sz = ht + botm if(iabs(ki).le.10.and.iabs(kj).le.10) then call mprism(s(i,j)/100., x, smxt, y, smyt, za, sz) else call mvline(s(i,j)/100., x, y, za, sz, am) endif f(ir+1,jr+1) = f(ir+1,jr+1) + calma(0.,0.,0.) 600 continue 500 continue 400 continue 300 continue call dpcent((j-1)*mxm+i, mxm*mym) 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)') * ixs, iys, mszx, mszy, mx, my, 9999.9, 0. do 700 j=1,my write(3,'((f7.1,9(1x,f7.1)))') (f(i,j),i=1,mx) 700 continue write(3,'(a8,i4,4x,4i8)') area, ncc, iorg,korg, ispa,ispb write(3,'(2i12,4i6,1x,f7.1,1x,f7.0)') * ixs, iys, mszx, mszy, mx, my, hnul, -1. do 800 j=1,my write(3,'((f7.1,9(1x,f7.1)))') (h(i,j),i=1,mx) 800 continue close(3) c if(flog(ldr:ldr).ne.' ') then open(99,file=flog,access='append') write(99,'(//a)') * 'CMAGF : calc. anomaly of fine model' write(99,'(a,a)') * ' Calc.Alt. data infile : ', fnam1(1:lrtrim(fnam1)) write(99,'(a,a)') * ' AMAG model infile : ', fnam2(1:lrtrim(fnam2)) write(99,'(a,a)') * ' Source Alt. infile : ', fnam3(1:lrtrim(fnam3)) write(99,'(a,a)') * ' Calc.Anom. outfile : ', fnam4(1:lrtrim(fnam4)) write(99,'(a,a,a)') '--------------- ', sdtm, '---------------' write(99,1000) area, ncc, iorg,korg, ispa,ispb, * xs, ys, mszx, mszy, mx, my, * xl, yl, mszxm,mszym, mxm,mym, * xst,yst, mszxt,mszyt, mxt,myt 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'/ * 3x, 'Calc. Alt. ', 2f10.2, 2i8, 2i6/ * 3x, 'AMAG model ', 2f10.2, 2i8, 2i6/ * 3x, 'Source Alt ', 2f10.2, 2i8, 2i6) 2000 format(' Half width of calc.window :', f6.2, * ' km (', i4, ' *', i4, ' calc.-mesh )') end