《重磁资料处理与解释》实验四位场边缘识别程序设计实验专业名称：地球物理学学生姓名：学生学号：指导老师：王万银、纪新林、纪晓琳、邱世灿提交日期：2016-1-31、基本原理地质目标体边缘时指断裂构造线、不同地质体的边界线，实际上是具有一定密度或磁性差异的地质体的边界线，在地质体的边缘附近，重、磁异常变化率较大，故所有的边缘识别方法均利用这一特点进行设计。

现在有重、磁位场边缘识别方法分为数理统计、数值计算以及其他三大类。

数值类边缘识别方法均利用极大值位置或零值位置确定地质体的边缘位置，其依据的理论基础是二度体铅垂台阶模型的重力异常特征。

在该模型边缘处重力异常总水平导数和解析信号振幅达到极大值、垂向导数达到零值。

故可以利用这些特征位置来确定二度体铅垂台阶的边缘位置，确定倾斜二度体、不规则二度体及三度体边缘位置的理论均为二度体铅垂台阶模型理论的推广，但确定的边缘位置和真实的位置有一定的偏差。

该偏差随着地质体边界形状、埋深、水平尺寸及物性差异等的变化而变化。

因此，边缘识别结果是一种定性或半定量解释结果，与定量解释结果有一定区别，识别结果可作为边缘位置定量反演的初值。

(1)垂向导数：垂向导数方法利用零值位置确定地质体的边缘位置，重力异常可以直接使用，对磁力异常必须转化为磁源重力异常或化极磁力异常才可以使用(,,)(,,)g x y zV D R x y zz∂=∂( 1.1)（2）解析信号振幅：解析信号振幅也是利用极大值位置来确定地质体的边缘位置适用于重、磁力异常(1.2)（3）总水平导数（THDR）(1.3)2、输入/输出数据格式设计依据上述原理，现在对上述各种边缘识别方法进行程序设计。

2.1 输入数据格式设计本次实验给了正演的重力异常数据，为.GRD格式，均为实型变量。

例如：DSAA201 201-1000.000000 1000.000000-1000.000000 1000.0000005.549671E-01 23.5398465.549671E-01 5.634658E-01 5.721339E-01 5.808522E-015.897312E-01 5.987253E-016.078691E-01 6.171604E-01…2.2 输出数据格式设计计算结果输出数据格式与输入格式对应，格式为.GRD格式，均为实型变量。

例如：DSAA201 201-1000.000 1000.000-1000.000 1000.000-0.1465084 0.3190881-3.6523044E-02 -3.3485338E-02 -3.3061244E-02 -3.2748722E-02 -3.2688729E-02-3.2654848E-02 -3.2723978E-02 -3.2787599E-02 -3.2927759E-02 -3.3138681E-02…2.3 参数文件数据格式设计将以上部分量保存在一个文件中，该文件名变量为cmdfile，字符串变量，长度不超过80，全路径名。

在该文件中保存的参数如下：输入数据文件名input_file，字符串变量，长度不超过80；输出vdr数据文件名output_file_vdr，字符串变量，长度不超过80；输出thdr数据文件名output_file_thdr，字符串变量，长度不超过80;输出asm数据文件名output_file_asm，字符串变量，长度不超过80factor_m：扩边比例因子，实型变量（>1）；3.总体设计此次程序采用IPO结构设计，首先通过读取cmd文件，得到相关输入参数：输入数据文件名gravity.grd、输出vdr文件名field_vrd.grd、输出thdr文件名field_thdr.grd、输出asm文件名field_asm.grd、扩边比例因子factor_m；然后确定确定扩边网格的大小，扩边数据点号位置；再从观测面位场数据文件中读取数据。

下一步，进行二维余弦扩边，将扩完边的数据进行快速二维傅里叶变换，转换到频率域；接下来在频率域求出在x,y,z方向的导数并反变换；最后求出VDR、THDR、ASM数据。

最后去除扩边部分后输出。

总体设计见表1。

4.测试结果分析：由图4.3可看出，VDR 方法的零值线可较准确识别模型体的边界；由图4.4可看出，ASM 的极大值点边界可大致识别模型体边界，但精度不是很高。

对比图4.2到4.5可以看出，THDR 方法对模型边界的识别效果是最好的。

5 结论及建议经测试，VDR 与THDR 对模型体的边界位置识别效果较好，而ASM 对模型体边界识别效果较差。

三种方法中，THDR 效果最好。

图4.3 垂向导数（VDR ）图4.2 重力异常原始图附录：边缘识别程序源代码!******************************************************************************************* *********! 程序功能：实现频率域位场导数运算进行边缘识别! cmd文件参数:! cmdfile:存放有关参数的文件名变量! input_file:观测面位场数据文件! output_file_vdr：场值的水平导数数据文件! output_file_thdr:场值垂向导数数据文件! output_file_asm:场值的解析信号振幅数据文件! factor_m：扩边因子! .grd文件参数：! N_point,N_line：点数(x方向)、线数(y方向)! x_min,x_max：x的最小值和最大值! y_min,y_max：y的最小值和最大值! Ur:初始观测面场值! 扩边参数：! m1,m2：x方向实际数据起点和终点点号位置! 1,m：x方向扩边后数据起点和终点点号位置! n1,n2：y方向实际数据起点和终点点号位置! 1,n3：y方向扩边后数据起点和终点点号位置! 求导参数：! field_re:初始观测面信号的实部! field_im：初始观测面信号的虚部! Px_re:x方向导数信号的实部! Px_im:x方向导数信号的虚部! Py_re:y方向导数信号的实部! Py_im:y方向导数信号的虚部! Pz_re:z方向导数信号的实部! Pz_im:z方向导数信号的虚部! W(m,n)：径向圆频率! field_vdr:对场值作水平导数的结果! field_thdr:对场值作垂向导数的结果! field_asm:场值的解析信号振幅的结果!******************************************************************************************* *********program deviationparameter(eigval=3.701411*1e5)character*(80)cmdfilecharacter*80 input_file,output_file_vdr,output_file_thdr,output_file_asmreal,allocatable:: field_re(:,:),field_im(:,:)real,allocatable:: Px_re(:,:),Px_im(:,:),Py_re(:,:),Py_im(:,:),Pz_re(:,:),Pz_im(:,:)real,allocatable:: field_vdr(:,:),field_thdr(:,:),field_asm(:,:)real,allocatable:: U(:),V(:),W(:,:)integer N_point,N_lineinteger m,n,m1,m2,n1,n2real factor_mreal xmin,xmax,ymin,ymax,dx,dycmdfile='deviation.cmd'call read_cmd(cmdfile,factor_m,input_file,output_file_vdr,output_file_thdr,output_file_asm)call read_grd(input_file,N_point,N_line,Xmin,Xmax,Ymin,Ymax)call calculate_mn(N_point,N_line,m,n,m1,m2,n1,n2,factor_m)allocate(field_re(1:m,1:n),field_im(1:m,1:n))allocate(Px_re(1:m,1:n),Px_im(1:m,1:n),Py_re(1:m,1:n),Py_im(1:m,1:n),Pz_re(1:m,1:n),Pz_im(1:m,1:n)) allocate(field_vdr(1:m,1:n),field_thdr(1:m,1:n),field_asm(1:m,1:n))allocate(U(1:m),V(1:n),W(1:m,1:n))call input_grd(field_re,input_file,m1,m2,n1,n2,m,n)call expand_cos_2D(m1,m2,m,n1,n2,n,field_re,field_im)call FFT2(field_re,field_im,m,n,2)CALL cal_dxdy(xmin,xmax,ymin,ymax,N_POINT,N_LINE,dx,dy)call WAVE2D(m,n,dx,dy,U,V,W)call factor(m,n,field_re,field_im,u,v,w,px_re,px_im,py_re,py_im,pz_re,pz_im)call FFT2(px_re,px_im,m,n,1)call FFT2(py_re,py_im,m,n,1)call FFT2(pz_re,pz_im,m,n,1)call deviration(m,n,px_re,py_re,pz_re,field_vdr,field_thdr,field_asm)call OUTPUT_GRD(field_vdr,output_file_vdr,m1,m2,n1,n2,m,n,eigval,xmin,xmax,ymin,ymax)call OUTPUT_GRD(field_thdr,output_file_thdr,m1,m2,n1,n2,m,n,eigval,xmin,xmax,ymin,ymax)call OUTPUT_GRD(field_asm,output_file_asm,m1,m2,n1,n2,m,n,eigval,xmin,xmax,ymin,ymax)deallocate(field_re,field_im,px_re,px_im,py_re,py_im,pz_re,pz_im,field_vdr,field_thdr,field_asm,u,v,w) end program!****************************************************************************! 子程序：read_cmd! 功能：读取参数文件! 输入参数说明：! cmdfile：参数文件名! 输出参数说明：! input_file:观测面位场数据文件! output_file_vdr：对场值作水平导数处理后的数据文件! output_file_thdr:对场值作垂向导数处理后的数据文件! output_file_asm:对场值作总导数处理后的数据文件! factor_m：扩边因子!**************************************************************************** Subroutine read_cmd(cmdfile,factor_m,input_file,output_file_vdr,output_file_thdr,output_file_asm) implicit nonecharacter*80 strcharacter*(*)cmdfilecharacter*(*)input_file,output_file_vdr,output_file_thdr,output_file_asmreal factor_mopen(10,file=cmdfile,status='old')read(10,*)str,input_fileread(10,*)str,output_file_vdrread(10,*)str,output_file_thdrread(10,*)str,output_file_asmread(10,*)str,factor_mclose(10)end Subroutine read_cmd!***************************************************************************! 子程序：read_grd! 功能：从原始观测.grd文件中读取相关参数! 输入参数说明：! filename_obser：输入文件名! 输出参数说明：! N_point,N_line：点数、线数! x_min,x_max：x的最小值和最大值! y_min,y_max：y的最小值和最大值!*************************************************************************** subroutine read_grd(input_file,N_point,N_line,Xmin,Xmax,Ymin,Ymax)implicit nonecharacter*(*)input_fileinteger N_point,N_linereal Xmin,Xmax,Ymin,Ymaxopen(10,file=input_file,status='old')Read(10,*)Read(10,*)N_line,N_pointRead(10,*)Xmin,XmaxRead(10,*)Ymin,YmaxClose(10)end subroutine read_grd!************************************************************************** ! 子程序：calculate_mn! 功能：确定扩边数据点号位置! 输入参数说明：! factor_m: 扩边比例因子（>1.0）! a,b：点数、线数! 输出参数说明：! m1,m2: x方向实际数据起点位置和终点位置点号! m:扩边后数据终点位置点号(起点位置点号为1)! n1,n2: y方向实际数据起点位置和终点位置点号! n:扩边后数据终点位置点号(起点位置点号为1)!**************************************************************************subroutine calculate_mn(a,b,m,n,m1,m2,n1,n2,factor_m)implicit noneinteger a,b,m,n,m1,m2,n1,n2integer mtemp,mu,nureal factor_mmtemp=aDO WHILE ((mod(mtemp,2).eq.0).and.(mtemp.ne.0))mtemp=mtemp/2End doIF (mtemp.eq.1) THENm=a*2ELSEmu=int(log(float(a))/0.693147+factor_m)m=2**muEND IFm1=1+(m-a)/2m2=m1+a-1write(*,*)m,apausemtemp=bDO WHILE ((mod(mtemp,2).eq.0).and.(mtemp.ne.0))mtemp=mtemp/2End doIF (mtemp.eq.1) THENn=b*2ELSEnu=int(log(float(b))/0.693147+factor_m)n=2**nuEND IFn1=1+(n-b)/2n2=n1+b-1write(*,*)m1,m2,n1,n2,m,npauseend subroutine calculate_mn!************************************************************************* ! 子程序：INPUT_GRD! 功能：读取grd文件中的数据! 输入参数说明：! filename_obser：输入文件名! m1,m2: x方向实际数据起点位置和终点位置点号! m:扩边后数据终点位置点号(起点位置点号为1)! n1,n2: y方向实际数据起点位置和终点位置点号! n:扩边后数据终点位置点号(起点位置点号为1)! 输出参数说明：! A：存放输出数据的二维数组名!************************************************************************* SUBROUTINE INPUT_GRD(A,input_file,m1,m2,n1,n2,m,n)character*(*)input_fileinteger m1,m2,n1,n2,m,nreal xmin,xmax,ymin,ymaxreal A(1:m,1:n)real i,j,kdo j=1,n,1do i=1,m,1A(i,j)=3.701411*1e10enddoenddoOpen(20,file=input_file,status='old')read(20,*)read(20,*)read(20,*)xmin,xmaxread(20,*)ymin,ymaxread(20,*)read(20,*) ((A(i,j),i=m1,m2),j=n1,n2)Close(20)END SUBROUTINE INPUT_GRD!************************************************************************* ! 子程序：expand_cos_2D! 功能：二维扩边子程序并为信号虚部赋值! 输入参数说明：! m1,m2: x方向实际数据起点位置和终点位置点号! m:扩边后数据终点位置点号(起点位置点号为1)! n1,n2: y方向实际数据起点位置和终点位置点号! n:扩边后数据终点位置点号(起点位置点号为1)! Ur:初始观测面信号的实部! Ui：初始观测面信号的虚部! 输出参数说明：! Ur:初始观测面信号的实部! Ui：初始观测面信号的虚部!*************************************************************************subroutine expand_cos_2D(m1,m2,m,n1,n2,n,Ur,Ui)implicit noneinteger m,n,m1,m2,n1,n2real Ur(1:m,1:n),Ui(1:m,1:n)real,allocatable::u(:),r(:)integer j,i,kallocate(u(1:m))do j=n1,n2,1do i=1,m,1u(i)=Ur(i,j)enddocall expand_cos_1d(1,m1,m2,m,u(1))do i=1,m,1Ur(i,j)=u(i)enddoenddodeallocate(u)allocate(r(1:n))do i=1,m,1do j=1,n,1r(j)=Ur(i,j)enddocall expand_cos_1d(1,n1,n2,n,r(1))do j=1,n,1Ur(i,j)=r(j)enddoenddodeallocate(r)do i=1,mdo j=1,nUi(i,j)=0enddoenddoend subroutine expand_cos_2D!************************************************************************** ! 子程序：expand_cos_1d! 功能：一维扩边子程序! 输入参数说明：! n0,n3：扩边后数据起点位置和终点位置! n1,n2：实际数据起点位置和终点位置! feild(i),(i=n1,n1+1,...,n2)：实际数据! 输出参数说明：! field(i),(i=n0,...,n1-1)：扩边后左边的数据! field(i),(i=n2+1,...,n3)：扩边后右边的数据!************************************************************************** Subroutine expand_cos_1d(n0,n1,n2,n3,Field)Real Field(n0:n3)pi=3.141592654Field (n0)=(Field (n1)+Field (n2))/2.0Field (n3)=Field (n0)i1=n0i2=n1DO i=i1,i2-1,1Field(i)=Field(i1)+cos(pi/2.0*(i2-i)/(i2-i1))*(Field(i2)-Field(i1)) End doi1=n2i2=n3DO i=i1+1,i2,1Field(i)=Field(i1)+cos(pi/2.0*(i2-i)/(i2-i1))*(Field(i2)-Field(i1)) End doEnd subroutine expand_cos_1d!************************************************************************ ! 功能：FFT2! 功能：复数组2-D快速Fourier变换! 输入参数说明：! m0,m3：x方向的起点和终点! n0,n3：y方向的起点和终点! field：输入信号(需要赋值给Freal，实部)! m,n：x,y方向扩边后数据终点点号位置（起始点号为1）! NF: 正、反变换标志量(1:反变换;2:正变换)! 输出参数说明：! Freal：信号的实部! Fimage：信号的虚部(对于实信号而言，赋值为0)! 对应频率分布说明：! 数据Freal(m,n)和Fimage(m,n)对应的频率分布位置为:! m 方向：0,1,.......,m/2-1,m/2,-(m/2-1),......,-1! n 方向：0,1,.......,n/2-1,n/2,-(n/2-1),......,-1!************************************************************************ SUBROUTINE FFT2(Freal,Fimage,m,n,nf)implicit noneINTEGER m,n,nfREAL Freal(1:m,1:n),Fimage(1:m,1:n)real,ALLOCATABLE::Treal(:),Timage(:)integer nmmax,ierr,i,jnmmax=max(m,n)allocate(Treal(1:nmmax),Timage(1:nmmax),STAT=ierr)if(ierr.ne.0) STOPDO i=1,m,1IF (n.ne.1) THENdo j=1,n,1Treal(j)=Freal(i,j)Timage(j)=Fimage(i,j)end docall FFT(Treal,Timage,n,nf)do j=1,n,1Freal(i,j)=Treal(j)Fimage(i,j)=Timage(j)end doEND IFEND DODO j=1,n,1IF(m.ne.1) THENdo i=1,m,1Treal(i)=Freal(i,j)Timage(i)=Fimage(i,j)end docall FFT(Treal,Timage,m,nf)do i=1,m,1Freal(i,j)=Treal(i)Fimage(i,j)=Timage(i)end doEND IFEND DOdeallocate(Treal,Timage,STAT=ierr)END SUBROUTINE FFT2!****************************************************************** ! 子程序：FFT! 功能：复数组1-D快速Fourier变换! 输入参数说明：! Xreal(n): 输入数据实部! Ximage(n): 输入数据虚部! N: 点数（N 必须为2 的幂次方）! NF: 正、反变换标志量(1:反变换;2:正变换)! 输出参数说明：! Xreal(n): 变换后频谱实部! Ximage(n): 变换后频谱虚部! 对应频率分布说明：! 数据Xreal(n)和Ximage(n)对应的频率分布位置为:! 0,1,.......,n/2-1,n/2,-(n/2-1),......,-1!***************************************************************** SUBROUTINE FFT(Xreal,Ximage,n,nf)implicit noneINTEGER n,nfREAL Xreal(1:n),Ximage(1:n)integer nu,n2,nu1,k,k1,k1n2,l,i,ibitrreal f,p,arg,c,s,treal,timagenu=int(log(float(n))/0.693147+0.001)n2=n/2nu1=nu-1f=float((-1)**nf)k=0DO l=1,nu,1DO while (k.lt.n)do i=1,n2,1p=ibitr(k/2**nu1,nu)arg=6.2831853*p*f/float(n)c=cos(arg)s=sin(arg)k1=k+1k1n2=k1+n2treal=Xreal(k1n2)*c+Ximage(k1n2)*stimage=Ximage(k1n2)*c-Xreal(k1n2)*sXreal(k1n2)=Xreal(k1)-trealXimage(k1n2)=Ximage(k1)-timageXreal(k1)=Xreal(k1)+trealXimage(k1)=Ximage(k1)+timagek=k+1end dok=k+n2END DOk=0nu1=nu1-1n2=n2/2END DODO k=1,n,1i=ibitr(k-1,nu)+1if(i.gt.k) thentreal=Xreal(k)timage=Ximage(k)Xreal(k)=Xreal(i)Ximage(k)=Ximage(i)Xreal(i)=trealXimage(i)=timageend ifEND DOIF(nf.ne.1) THENdo i=1,n,1Xreal(i)=Xreal(i)/float(n)Ximage(i)=Ximage(i)/float(n)end doEND IFEND SUBROUTINE FFTFUNCTION IBITR(J,NU)implicit noneinteger ibitr,j,nuinteger j1,itt,i,j2j1=jitt=0do i=1,nu,1j2=j1/2itt=itt*2+(j1-2*j2)ibitr=ittj1=j2end doEND FUNCTION IBITR!*************************************************************************** ! 子程序：cal_dxdy! 功能：计算x,y方向的点距! 输入参数说明：! x_min,x_max：x的最小值和最大值! y_min,y_max：y的最小值和最大值! N_point,N_line：点数(x方向)、线数(y方向)! 输出参数说明：! dx,dy：x,y方向的点距!*************************************************************************** subroutine cal_dxdy(xmin,xmax,ymin,ymax,N_POINT,N_LINE,dx,dy)implicit nonereal xmin,xmax,ymin,ymaxinteger N_POINT,N_LINEreal dx,dydx=(xmax-xmin)/N_POINTdy=(ymax-ymin)/N_LINEend subroutine cal_dxdy!******************************************************************! 子程序：WAVE2D! 功能：计算2D径向圆频率W! 输入参数说明：! dx: x 方向点距! dy: y 方向线距! m: 点数（M 必须为2 的幂次方）! n: 线数（N 必须为2 的幂次方）! 输出参数说明：! W(m,n): 径向圆频率!******************************************************************SUBROUTINE WAVE2D(m,n,dx,dy,U,V,W)implicit noneINTEGER m,nREAL dx,dyREAL W(1:m,1:n),U(1:m),V(1:n)real pi,delx,delyinteger midm,midn,i,j,xx,yymidm=m/2+1midn=n/2+1delx=float(m)/dxdely=float(n)/dydo j=1,n,1yy=jif(yy.gt.midn) yy=yy-nv(j)=dely*(yy-1)end dodo i=1,m,1xx=iif(xx.gt.midm) xx=xx-mu(i)=delx*(xx-1)end dodo j=1,n,1do i=1,m,1w(i,j)=sqrt(u(i)*u(i)+v(j)*v(j))end doend doEND SUBROUTINE WAVE2D!*********************************************************************************** ! 子程序：factor! 功能：计算x,y,z方向导数的实部和虚部!! 输入参数说明：! m: 点数（M 必须为2 的幂次方）! n: 线数（N 必须为2 的幂次方）! field_re:初始观测面信号的实部! field_im：初始观测面信号的虚部! W(m,n)：径向圆频率! 输出参数说明：! Px_re:x方向导数信号的实部! Px_im:x方向导数信号的虚部! Py_re:y方向导数信号的实部! Py_im:y方向导数信号的虚部! Pz_re:z方向导数信号的实部! Pz_im:z方向导数信号的虚部!*********************************************************************************** subroutine factor(m,n,field_re,field_im,u,v,w,px_re,px_im,py_re,py_im,pz_re,pz_im)implicit noneinteger m,nreal field_re(1:m,1:n),field_im(1:m,1:n)real px_re(1:m,1:n),px_im(1:m,1:n),py_re(1:m,1:n),py_im(1:m,1:n),pz_re(1:m,1:n),pz_im(1:m,1:n) real U(1:m),V(1:n),W(1:m,1:n)integer i,jreal pipi=3.1415926do i=1,m,1do j=1,n,1px_re(i,j)=field_im(i,j)*u(i)*(-1)*pi*2.0px_im(i,j)=field_re(i,j)*u(i)*pi*2.0py_re(i,j)=field_im(i,j)*v(j)*(-1)*pi*2.0py_im(i,j)=field_re(i,j)*v(j)*pi*2.0pz_re(i,j)=field_re(i,j)*w(i,j)*pi*2.0pz_im(i,j)=field_im(i,j)*w(i,j)*pi*2.0end doend doend subroutine factor!***************************************************************************! 子程序：deviration! 功能：计算异常的水平导数、垂向导数及总导数!! 输入参数说明：! m: 点数（M 必须为2 的幂次方）! n: 线数（N 必须为2 的幂次方）! Px_re:x方向导数信号的实部! Px_im:x方向导数信号的虚部! Py_re:y方向导数信号的实部! Py_im:y方向导数信号的虚部! Pz_re:z方向导数信号的实部! Pz_im:z方向导数信号的虚部! 输出参数说明：! field_vdr:对场值作水平导数的结果! field_thdr:对场值作垂向导数的结果! field_asm:对场值作总导数的结果!************************************************************************** subroutine deviration(m,n,px_re,py_re,pz_re,field_vdr,field_thdr,field_asm)implicit noneinteger m,nreal px_re(1:m,1:n),py_re(1:m,1:n),pz_re(1:m,1:n)real field_vdr(1:m,1:n),field_thdr(1:m,1:n),field_asm(1:m,1:n)integer i,jdo i=1,m,1do j=1,n,1field_vdr(i,j)=pz_re(i,j)field_thdr(i,j)=sqrt(px_re(i,j)**2+py_re(i,j)**2)field_asm(i,j)=sqrt(field_vdr(i,j)**2+field_thdr(i,j)**2)end doend doend subroutine deviration!*************************************************************************** ! 子程序：OUTPUT_GRD! 功能：输出结果!! 输入参数说明：! A：存放输出数据的二维数组! m: 点数（M 必须为2 的幂次方）! n: 线数（N 必须为2 的幂次方）! m1,m2：实际数据x方向起点位置和终点位置! n1,n2：实际数据y方向起点位置和终点位置! 输出参数说明：! x_min,x_max：x的最小值和最大值! y_min,y_max：y的最小值和最大值! filename:存放输出结果的文件名! field_vdr:对场值作水平导数的结果! field_thdr:对场值作垂向导数的结果! field_asm:场值的解析信号振幅的结果!************************************************************************** SUBROUTINE OUTPUT_GRD(A,filename,m1,m2,n1,n2,m,n,eigval,xmin,xmax,ymin,ymax) real eigvalinteger m1,m2,n1,n2,m,nreal A(1:m,1:n)character*(*) filenamereal amin,amaxamin=HUGE(amin)amax=-HUGE(amax)write(*,*)m1,m2,n1,n2,m,npauseeigval=3.701411*1e5DO j=n1,n2,1do i=m1,m2,1if(ABS(A(i,j)).lt.eigval) thenamin=MIN(amin,A(I,j))amax=MAX(amax,A(I,j))end ifend doEND DOOPEN(20,file=filename,status='unknown',form='formatted') write (20,'(A)')'DSAA'write (20,*)m2-m1+1,n2-n1+1write (20,*)xmin,xmaxwrite (20,*)ymin,ymaxwrite(20,*)amin,amaxDo j=n1,n2,1write(20,*) (A(i,j),i=m1,m2)End doClose(20)END subroutine OUTPUT_GRD。