c--------( AMAGC : inversion until converge )-------- parameter (LHW=100, LXO=1201, LYO=1201) parameter (LXM=LXO+LHW+LHW, LYM=LYO+LHW+LHW) c dimension g((LHW*2+1)*(LHW*2+1)), coef(3) dimension f(LXO,LYO), r(LXO,LYO), q(LXO,LYO) dimension s(LXM,LYM), p(LXM,LYM), v(LXM,LYM) character*80 wdr, flog, fnam1, fnam2, fnam3, fnam4 character buf*80, area*8, area1*8, out*72, salt*8, yn*1, sdtm*20 c call premsg('AMAGC : inversion until converge') 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)) call gparma('Enter Mag.Anomaly data filename ==> ', * 81-ldr, fnam1(ldr:80)) open(12,file=fnam1,status='old') 1 read(12,'(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(12,'(a)') buf read(buf,*) ixs, iys, mszx, mszy, mx, my, vnul, alt 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(12,*) ((f(i,j),i=1,mx),j=1,my) close(12) call prompt(' Remove Linear trend ? ') call gparma(' "y" (Yes) or "n" (No) ==> ', 1, yn) if(yn.ne.'y'.and.yn.ne.'n') call abendm('invalid parameter') c call gparma('Enter COEF matrix data filename ==> ', * 81-ldr, fnam2(ldr:80)) open(20,file=fnam2,form='unformatted',status='old') read(20) area1, ncc1, iorg1,korg1, ispa1,ispb1, * ixs1, iys1, mszx1, mszy1, mx1, my1, alt1, * ixl, iyl, mszxm, mszym, mxm, mym, nhx, nhy if(area1.ne.area.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(ixs1.ne.ixs.or.iys1.ne.iys .or. * mszx1.ne.mszx.or.mszy1.ne.mszy .or. * mx1.ne.mx.or.my1.ne.my) * call abendm('Header inconsistent') if(alt.lt.0.) alt = alt1 write(salt,'(1x,f7.0)') alt if (alt.eq.0.) salt = ' var.' if (alt.lt.0.) salt = ' undef' xs = float(ixs)/1000. ys = float(iys)/1000. xl = float(ixl)/1000. yl = float(iyl)/1000. 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,a8)') * '--- mag. obs. ', xs, ys, mszx, mszy, mx, my, salt call premsg(out) write(out,'(a,2f10.2,2i8,2i6,a8)') * '--- source model', xl,yl, mszxm,mszym, mxm,mym, ' --- ' call premsg(out) if(mxm.gt.LXM.or.mym.gt.LYM) call abendm('too large array size') hw = float(mszx*nhx) / 1000. nfx = nhx + 1 + nhx nfy = nhy + 1 + nhy nnfw = nfx * nfy c call gparma('Enter AMAG model in/out filename ==> ', * 81-ldr, fnam3(ldr:80)) open(21,file=fnam3,status='old') 2 read(21,'(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(21,*) ixl1, iyl1, mszxm1, mszym1, mxm1, mym1, anul, dmy, klp if(area1.ne.area.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.ixl.or.iyl1.ne.iyl .or. * mszxm1.ne.mszxm.or.mszym1.ne.mszym .or. * mxm1.ne.mxm.or.mym1.ne.mym) * call abendm('Header inconsistent') read(21,*) ((s(i,j),i=1,mxm),j=1,mym) rewind(21) c call premsg('Enter Auxiliary(2) output filename') call gparma(' ( blank to suppress output ) ==> ', * 81-ldr, fnam4(ldr:80)) kout = 0 if(fnam4(ldr:ldr).ne.' ') then open(8,file=fnam4,status='new') kout = 1 endif call premsg('Threshold percentage to regard converged') call gparmf(' (cont. 5 loops) [default: 2%] ==> ', thrsh) if (thrsh.lt.0. .or. thrsh.ge.20.) * call abendm('invalid parameter') if (thrsh.eq.0.) thrsh = 2. call clspin() call strdtm(sdtm) c call dpcini('.... read COEF matrix (0) : ') do 100 j=1,mym do 100 i=1,mxm read(20) ip, jp, k, k, (c,k=1,nnfw) if(ip.ne.i.or.jp.ne.j) call abendm('COEF file broken') p(i,j) = 0. call dpcent((j-1)*mxm+i, mxm*mym) 100 continue rewind(20) c write(6,'(a)') '----------------------------------------' write(6,'(a)') 'AMAGC : inversion until converge' write(6,'(a,a)') * ' Mag. Anom. infile : ', fnam1(1:lrtrim(fnam1)) write(6,'(a,a)') * ' COEF matrix infile : ', fnam2(1:lrtrim(fnam2)) write(6,'(a,a)') * ' AMAG model i/o file : ', fnam3(1:lrtrim(fnam3)) if(kout.eq.1) write(6,'(a,a)') * ' Aux(2) output file : ', fnam4(1:lrtrim(fnam4)) write(6,1000) area, ncc, iorg,korg, ispa,ispb, * xs, ys, mszx, mszy, mx, my, salt, * xl, yl, mszxm, mszym, mxm, mym, ' --- ' write(6,2000) hw, nhx, nhy c n = 0 ss = 0. do 200 j=1,my do 200 i=1,mx if(f(i,j).eq.vnul) goto 200 ss = ss + f(i,j) n = n + 1 200 continue ft = ss / float(n) ss = 0. do 220 j=1,mym do 220 i=1,mxm ss = ss + s(i,j)*s(i,j) 220 continue ds = sqrt(ss / float(mxm*mym)) c if (yn.eq.'y') then cmx = float(mx+1) / 2. cmy = float(my+1) / 2. call sm2opn(1, 1) do 400 j=1,my y = float(j) - cmy do 400 i=1,mx x = float(i) - cmx if(f(i,j).ne.vnul) call sm2ex(x, y, f(i,j)-ft) 400 continue call sm2cls(coef, 3, 1) cf = coef(1) + ft af = coef(2) bf = coef(3) write(6,3000) cf, af, bf do 500 j=1,my y = float(j) - cmy do 500 i=1,mx x = float(i) - cmx if(f(i,j).ne.vnul) f(i,j) = f(i,j) - (af*x + bf*y + cf) 500 continue if(kout.eq.1) then write(8,'(a8,i4,4x,4i8)') * 'Obs-Trnd', ncc, iorg,korg, ispa,ispb write(8,'(2i12,4i6,1x,f7.1,1x,f7.0)') * ixs,iys, mszx,mszy, mx,my, vnul, -1. do 550 j=1,my write(8,'((f7.1,9(1x,f7.1)))') (f(i,j),i=1,mx) 550 continue write(6,'(a)') * 'AUX2 <== Obs-Trnd : Trend-removed Observation' endif else af = 0. bf = 0. cf = 0. endif c call conv1(g, nfx, nfy, LXM, s, mxm, mym, LXO, r, mx, my) n = 0 ss = 0. do 10 j=1,my do 10 i=1,mx if(f(i,j).eq.vnul) then r(i,j) = 0. else r(i,j) = f(i,j) - r(i,j) ss = ss + r(i,j)*r(i,j) n = n + 1 endif 10 continue dv = sqrt(ss / float(n)) lp = 0 write(6,4000) klp, lp, ds, dv open(9,status='scratch') write(9,*) klp, lp, ds, dv c pvv = 1. 7 klp = klp + 1 rv = dv pvvp = pvv call conv2(g, nfx, nfy, LXO, r, mx, my, LXM, v, mxm, mym) pvv = 0. do 20 j=1,mym do 20 i=1,mxm pvv = pvv + v(i,j)*v(i,j) 20 continue a = pvv / pvvp do 30 j=1,mym do 30 i=1,mxm p(i,j) = v(i,j) + a*p(i,j) 30 continue call conv1(g, nfx, nfy, LXM, p, mxm, mym, LXO, q, mx, my) pqq = 0. do 40 j=1,my do 40 i=1,mx pqq = pqq + q(i,j)*q(i,j) 40 continue a = pvv / pqq ss = 0. do 50 j=1,mym do 50 i=1,mxm s(i,j) = s(i,j) + a*p(i,j) ss = ss + p(i,j)*p(i,j) 50 continue ds = a * sqrt(ss / float(mxm*mym)) n = 0 ss = 0. do 60 j=1,my do 60 i=1,mx if(f(i,j).ne.vnul) then r(i,j) = r(i,j) - a*q(i,j) ss = ss + r(i,j)*r(i,j) n = n + 1 endif 60 continue dv = sqrt(ss / float(n)) if (((rv-dv)/rv)*100. .ge. thrsh) then lp = 0 else lp = lp + 1 endif write(6,4000) klp, lp, ds, dv write(9,*) klp, lp, ds, dv if (lp.lt.5.and.dv.ge.0.1) goto 7 c write(6,'(1x)') endfile(9) rewind(9) c write(21,'(a8,i4,4x,4i8)') area, ncc, iorg,korg, ispa,ispb write(21,'(2i12,4i6,1x,f7.1,1x,f7.0,i6)') * ixl, iyl, mszxm, mszym, mxm, mym, anul, -1., klp do 700 j=1,mym do 710 i=1,mxm if(s(i,j).lt.-9999.8) s(i,j) = -9999.8 if(s(i,j).gt.+9999.8) s(i,j) = +9999.8 710 continue write(21,'((f7.1,9(1x,f7.1)))') (s(i,j),i=1,mxm) 700 continue close(21) if(kout.eq.1) then call conv1(g, nfx, nfy, LXM, s, mxm, mym, LXO, r, mx, my) do 800 j=1,my y = float(j) - cmy do 800 i=1,mx x = float(i) - cmx if(f(i,j).eq.vnul) then r(i,j) = vnul q(i,j) = vnul else q(i,j) = r(i,j) + (af*x + bf*y + cf) f(i,j) = f(i,j) - r(i,j) endif 800 continue write(8,'(a8,i4,4x,4i8)') 'SynthMdl', ncc, iorg,korg, ispa,ispb write(8,'(2i12,4i6,1x,f7.1,1x,f7.0)') * ixs,iys, mszx,mszy, mx,my, vnul, -1. do 810 j=1,my write(8,'((f7.1,9(1x,f7.1)))') (r(i,j),i=1,mx) 810 continue write(6,'(a)') 'AUX2 <== SynthMdl : Synthetic model anomaly ' if (yn.eq.'y') then write(8,'(a8,i4,4x,4i8)') * 'Mdl+Trnd', ncc, iorg,korg, ispa,ispb write(8,'(2i12,4i6,1x,f7.1,1x,f7.0)') * ixs,iys, mszx,mszy, mx,my, vnul, -1. do 820 j=1,my write(8,'((f7.1,9(1x,f7.1)))') (q(i,j),i=1,mx) 820 continue write(6,'(a)') 'AUX2 <== Mdl+Trnd : Trend-added Synth. model' endif write(8,'(a8,i4,4x,4i8)') 'Residual', ncc, iorg,korg, ispa,ispb write(8,'(2i12,4i6,1x,f7.1,1x,f7.0)') * ixs,iys, mszx,mszy, mx,my, vnul, -1. do 830 j=1,my write(8,'((f7.1,9(1x,f7.1)))') (f(i,j),i=1,mx) 830 continue write(6,'(a)') * 'AUX2 <== Residual : Residual <(Obs.) - (Model)>' close(8) endif close(20) c if(flog(ldr:ldr).ne.' ') then open(99,file=flog,access='append') write(99,'(//a)') 'AMAGC : inversion until converge' write(99,'(a,a)') * ' Mag. Anom. infile : ', fnam1(1:lrtrim(fnam1)) write(99,'(a,a)') * ' COEF matrix infile : ', fnam2(1:lrtrim(fnam2)) write(99,'(a,a)') * ' AMAG model i/o file : ', fnam3(1:lrtrim(fnam3)) if(kout.eq.1) write(99,'(a,a)') * ' Aux(2) output file : ', 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, salt, * xl, yl, mszxm, mszym, mxm, mym, ' --- ' write(99,2000) hw, nhx, nhy write(99,3000) cf, af, bf 92 read(9,*,end=93) klp, lp, ds, dv write(99,4000) klp, lp, ds, dv goto 92 93 call strdtm(sdtm) write(99,'(a,a,a)') '=============== ', sdtm, '===============' close(99) endif close(9) stop c 1000 format('Areaname : ', a8, 5x, 'Coord.', i4, 4x, 4i8/ * 3x, ' (mesh) ', ' xs ', ' ys ', * ' msz-x', ' msz-y', ' mx', ' my', ' alt '/ * 3x, 'mag. obs. ', 2f10.2, 2i8, 2i6, a8/ * 3x, 'source model', 2f10.2, 2i8, 2i6, a8) 2000 format(' Half width of calc.window :', f6.2, ' km (', i4, * ' *', i4, ' mag.mesh )') 3000 format(' Trend surface in nT :', f8.2, ' + (', f7.3, ' * x) + (', * f7.3, ' * y)'/ 5x, '( x/y : N/E coord. in mesh unit ', * 'referred to the center of area )') 4000 format(' loop', i4, I3, ' rms.Step :', f9.3, * ' rms.Residual :', f7.2) end c--- subroutine conv1(g, nfx, nfy, lsx, s, mxm, mym, lrx, r, mx, my) dimension g(nfx,nfy) dimension s(lsx,*), r(lrx,*) character*8 area nhx = (nfx-1) / 2 nhy = (nfy-1) / 2 call dpcini('.... read COEF matrix (1) : ') read(20) area, (i,k=1,20) do 1 jp=1,my do 1 ip=1,mx r(ip,jp) = 0. 1 continue do 2 js=1,mym do 2 is=1,mxm read(20) k, k, kic, kjc, g kis = kic - nhx kjs = kjc - nhy do 3 j=1,nfy jp = kjs + j if(jp.le.0.or.jp.gt.my) goto 3 do 4 i=1,nfx ip = kis + i if(ip.le.0.or.ip.gt.mx) goto 4 r(ip,jp) = r(ip,jp) + g(i,j)*s(is,js) 4 continue 3 continue call dpcent((js-1)*mxm+is, mxm*mym) 2 continue rewind(20) return end c--- subroutine conv2(g, nfx, nfy, lrx, r, mx, my, lsx, s, mxm, mym) dimension g(nfx,nfy) dimension s(lsx,*), r(lrx,*) character*8 area nhx = (nfx-1) / 2 nhy = (nfy-1) / 2 call dpcini('.... read COEF matrix (2) : ') read(20) area, (i,k=1,20) do 1 js=1,mym do 1 is=1,mxm read(20) k, k, kic, kjc, g kis = kic - nhx kjs = kjc - nhy ss = 0. do 2 j=1,nfy jp = kjs + j if(jp.le.0.or.jp.gt.my) goto 2 do 3 i=1,nfx ip = kis + i if(ip.le.0.or.ip.gt.mx) goto 3 ss = ss + g(i,j)*r(ip,jp) 3 continue 2 continue s(is,js) = ss call dpcent((js-1)*mxm+is, mxm*mym) 1 continue rewind(20) return end