c--------( CALMAS : calc. model anomaly on a given surface )-------- parameter (LXO=1201, LYO=1201, MNM=20) parameter (DPI=3.14159*2.) c dimension f(LXO,LYO), h(LXO,LYO) dimension pv(7,MNM), pam(MNM), pth(MNM), pph(MNM), ipv(MNM) character*80 wdr, flog, fnam1, fnam2 character buf*80, out*72, area*8, area1*8, area2*8, sdtm*20 character*5 nmdl(4)/ 'Block', 'hRect', 'vLine', 'Point'/ c call premsg('CALMAS : calc. model anomaly on a given surface') call opnpin() ldr = lwkdir(50,wdr) + 1 flog = wdr fnam1 = wdr fnam2 = wdr call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80)) call gparma('Enter Calc.Surface 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. do 100 j=1,my do 100 i=1,mx f(i,j) = 0. 100 continue ixh = ixl + mszx*(mx-1) iyh = iyl + mszy*(my-1) call gparma('Enter Calc. MA output filename ==> ', * 81-ldr, fnam2(ldr:80)) open(2,file=fnam2,status='new') c write(6,'(a)') '----------------------------------------' write(6,'(a)') 'CALMAS : calc. model anomaly on a given surface' write(6,'(a,a)') * ' Calc. alt. data infile : ', fnam1(1:lrtrim(fnam1)) write(6,'(a,a)') * ' Model anomaly outfile : ', fnam2(1:lrtrim(fnam2)) write(6,1000) area, ncc, iorg,korg, ispa,ispb, * xl, yl, mszx, mszy, mx, my, * xl, yl, mszx, mszy, mx, my, ' var.' c call premsg('') call premsg('How to specify the source location') call premsg(' 0: Coordinates values or') call gparmi(' 1: Distance from SW corner ==> ', kl) if(kl.lt.0.or.kl.gt.1) call abendm('invalid selection') call premsg('Specify Ambient field Direction') call gparmf2(' Inclin., Declin. (in degrees) : ', dip, dec) if(dip.lt.-90..or.dip.gt.90. .or. * dec.lt.-180..or.dec.ge.360.) call abendm('invalid value') c nm = 0 2 call premsg('Select source model type:') call prompt(' 1:Block 2:hRect 3:vLine 4:Point or 0:End') call gparmi(' ==> ', ks) if(ks.lt.0.or.ks.gt.4) call abendm('invalid selection') if(ks.eq.0) goto 3 nm = nm + 1 if(nm.gt.MNM) call abendm('too many models') ipv(nm) = ks call premsg(' Horizontal location:') if(ks.eq.1.or.ks.eq.2) then call prompt(' Northing Center (cx), N-S size (rx) (in km)') call gparmf2(' ==> ', cx, rx) if(kl.eq.1) cx = cx + xl pv(1,nm) = cx pv(2,nm) = rx call prompt(' Easting Center (cy), E-W size (ry) (in km) ') call gparmf2(' ==> ', cy, ry) if(kl.eq.1) cy = cy + yl pv(3,nm) = cy pv(4,nm) = ry else call gparmf(' Northing (xc) (in km) ==> ', xc) if(kl.eq.1) xc = xc + xl pv(1,nm) = xc pv(2,nm) = 0. call gparmf(' Easting (yc) (in km) ==> ', yc) if(kl.eq.1) yc = yc + yl pv(3,nm) = yc pv(4,nm) = 0. endif call premsg(' Depth below Sea-level:') if(ks.eq.1.or.ks.eq.3) then call prompt(' Top depth (za), and Thickness (rz) (in km)') call gparmf2(' ==> ', za, rz) pv(5,nm) = za pv(6,nm) = rz else call gparmf(' (in km) ==> ', zc) pv(5,nm) = zc pv(6,nm) = 0. endif if(ks.eq.4) then call prompt(' Equivalent source volume (in cubic-km)') call gparmf(' ==> ', vol) else if(ks.eq.3) then call prompt(' Equivalent source area (in square-km)') call gparmf(' ==> ', vol) else if(ks.eq.2) then call prompt(' Equivalent source thickness (in km)') call gparmf(' ==> ', vol) else vol = 0. endif pv(7,nm) = vol call gparmf(' Magnetization Intensity (A/m) : ', am) pam(nm) = am call premsg(' Magnetization Direction ?') call gparmf2(' Inclin., Declin. (in degrees) : ', the, phi) if(the.lt.-90..or.the.gt.90. .or. * phi.lt.-180..or.phi.ge.360.) call abendm('invalid value') pth(nm) = the pph(nm) = phi write(6,2000) nm, nmdl(ipv(nm)), pam(nm), pth(nm), pph(nm), * (pv(i,nm),i=1,7) c call magafd(the, phi, dip, dec) if(ks.eq.4) then call mpoint(am, xc, yc, zc, vol) else if(ks.eq.3) then call mvline(am, xc, yc, za, rz, vol) else if(ks.eq.2) then call mhrect(am, cx, rx, cy, ry, zc, vol) else call mprism(am, cx, rx, cy, ry, za, rz) endif 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 + mszx*(i-1)) / 1000. yp = float(iyl + mszy*(j-1)) / 1000. zp = - h(i,j) / 1000. f(i,j) = f(i,j) + calma(xp, yp, zp) endif 200 continue goto 2 c 3 call clspin() write(2,'(a8,i4,4x,4i8)') area, ncc, iorg,korg, ispa,ispb write(2,'(2i12,4i6,1x,f7.1,1x,f7.0)') * ixl, iyl, mszx, mszy, mx, my, 9999.9, 0. do 300 j=1,my write(2,'((f7.1,9(1x,f7.1)))') (f(i,j),i=1,mx) 300 continue write(2,'(a8,i4,4x,4i8)') area, ncc, iorg,korg, ispa,ispb write(2,'(2i12,4i6,1x,f7.1,1x,f7.0)') * ixl, iyl, mszx, mszy, mx, my, hnul, -1. do 400 j=1,my write(2,'((f7.1,9(1x,f7.1)))') (h(i,j),i=1,mx) 400 continue close(2) c if(flog(ldr:ldr).ne.' ') then open(99,file=flog,access='append') write(99,'(//a)') * 'CALMAS : calc. model anomaly on a given surface' write(99,'(a,a)') * ' Calc. alt. data infile : ', fnam1(1:lrtrim(fnam1)) write(99,'(a,a)') * ' Model anomaly outfile : ', fnam2(1:lrtrim(fnam2)) write(99,1000) area, ncc, iorg,korg, ispa,ispb, * xl, yl, mszx, mszy, mx, my, * xl, yl, mszx, mszy, mx, my, ' var.' do 500 n=1,nm write(99,2000) n, nmdl(ipv(n)), pam(n), pth(n), pph(n), * (pv(i,n),i=1,7) 500 continue 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, 'calc. surf. ', 2f10.2, 2i8, 2i6/ * 3x, 'model anom. ', 2f10.2, 2i8, 2i6, a8) 2000 format(i4, 2x, a5, f8.2, 2f8.1/ 6x, 3(f9.2, f7.2), f9.4) end