subroutine mlayer(mag,nx,ny,dx,dy,z1,z2,mi,md,fi,fd,store) c c Subroutine MLAYER calculates the total-field anomaly on a two- c dimensional grid due to a horizontal layer with two- c dimensional magnetization. The following steps are involved: c (1) Fourier transform the magnetization, (2) multiply by the c earth filter, and (3) inverse Fourier transform the product. c Magnetization is specified on a rectangular grid with x and y c axes directed north and east, respectively. Z axis is down. c Distance units irrelevant but must be consistent. Requires c subroutines FOURN, DIRCOS, KVALUE, and MFILT. c c Input parameters: c nx - number of elements in the sout_to_north direction. c ny - number of elements in the west-to-east direction. c (NOTE: both nx and ny must be a power of two.) c mag - a singly dimensioned real array containing the c two-dimensional magnetization (in A/m). Elements should c be in order of west to east, then south to north (i.e., c element 1 is southwest corner, element ny is c southeast corner, element (nx-1)*ny+1 is northwest c corner, and element ny*nx is northeast corner. c store - a singly dimensioned real array used internally. c It should be dimensioned at least 2*nx*ny. c dx - sample interval in the x direction. c dy - sample interval in the y direction. c z1 - depth to top of layer. Must be > 0. c z2 - depth to bottom of layer. Must be > z1. c mi - inclination of magnetization, in degrees positive below c horizontal. c md - declination of magnetization, in degrees east of north. c fi - inclination of regional field. c fd - declination of regional field. c c Output parameters: c mag - upon output, mag contains the total-field anomaly c (in nT) with same orientation as above. c complex cmag,mfilt,cmplx real mag,mi,md,mx,my,mz,kx,ky dimension mag(nx*ny),store(2*nx*ny),nn(2) data pi/3.14159265/,t2nt/1.e9/ index(i,j,ncol)=(j-1)*ncol+i nn(1)=ny nn(2)=nx ndim=2 call dircos(mi,md,0.,mx,my,mz) call dircos(fi,fd,0.,fx,fy,fz) dkx=2.*pi/(nx*dx) dky=2.*pi/(ny*dy) do 10 j=1,nx do 10 i=1,ny ij=index(i,j,ny) store(2*ij-1)=mag(ij) 10 store(2*ij)=0. call fourn(store,nn,ndim,-1) do 20 j=1,nx do 20 i=1,ny ij=index(i,j,ny) call kvalue(i,j,nx,ny,dkx,dky,kx,ky) cmag=cmplx(store(2*ij-1),store(2*ij)) cmag=cmag*mfilt(kx,ky,mx,my,mz,fx,fy,fz,z1,z2) store(2*ij-1)=real(cmag) 20 store(2*ij)=aimag(cmag) call fourn(store,nn,ndim,+1) do 30 j=1,nx do 30 i=1,ny ij=index(i,j,ny) 30 mag(ij)=store(2*ij-1)*t2nt/(nx*ny) return end