重磁异常空间域正演.docx
《重磁异常空间域正演.docx》由会员分享，可在线阅读，更多相关《重磁异常空间域正演.docx（17页珍藏版）》请在冰豆网上搜索。
重磁异常空间域正演
!
对给定的密度体和磁性体计算起伏观测面以及平面上的重力异常、磁力异常以及化极磁力异常。
PROGRAMComplex_body_forward
IMPLICITNONE
CHARACTER*80cmdfile,sourcefile,stationfile,output_grafile,
&output_magfile,output_polfile
INTEGERcomponent,order,plane!
场的分量标识，导数的阶数，观测面的形态标识
INTEGERn_source,n_station,i,j,mpoint,line,cs,cxyz
REALxmin,xmax,ymin,ymax
REAL,ALLOCATABLE:
:
SOURCE（:
:
）,XYZ（:
:
）
REAL,ALLOCATABLE:
:
ANOM_GRA（:
）,ANOM_MAG（:
）,ANOM_POL（:
）
cmdfile='cmd.dat'
cs=11
cxyz=3
CALLINPUT_CMDFILE（component,order,plane,cmdfile,sourcefile,
&stationfile,output_grafile,output_magfile,output_polfile）
CALLCHECK_CMD（component,order,plane,stationfile）
CALLGET_NUMBER_SOURCE（sourcefile,n_source）
ALLOCATE（SOURCE（1:
n_source,1:
cs））
CALLINPUT_SOURCE（sourcefile,SOURCE,n_source,cs）
CALLGET_n_station（stationfile,plane,xmin,xmax,ymin,ymax,
&mpoint,line,n_station）
ALLOCATE（XYZ（1:
n_station,1:
cxyz））
CALLINPUT_STATION（stationfile,plane,xmin,xmax,ymin,ymax,
&mpoint,line,XYZ,n_station）
ALLOCATE（ANOM_GRA（1:
n_station）,ANOM_MAG（1:
n_station））
ALLOCATE（ANOM_POL（1:
n_station））
CALLCAL_ANOM（SOURCE,n_source,cs,XYZ,n_station,cxyz,
&ANOM_GRA,ANOM_MAG,ANOM_POL,component）
CALLOUTPUT（plane,XYZ,ANOM_GRA,ANOM_MAG,ANOM_POL,n_station,xmin
&,xmax,ymin,ymax,mpoint,line,output_grafile,output_magfile
&,output_polfile）
ENDPROGRAM
!
读取场的分量标识、导数的阶数、观测面的形态标识
!
读取场源参数文件名、测点文件名、输出文件名
SUBROUTINEINPUT_CMDFILE（component,order,plane,cmdfile,sourcefile,
&stationfile,output_grafile,output_magfile,output_polfile）
IMPLICITNONE
INTEGERcomponent,order,plane
CHARACTER*80cmdfile,sourcefile,stationfile,output_grafile,
&output_magfile,output_polfile
OPEN（12,FILE=cmdfile,ACTION='READ'）
READ（12,*）output_polfile
READ（12,*）output_polfile
READ（12,*）output_polfile
READ（12,*）component
READ（12,*）output_polfile
READ（12,*）order
READ（12,*）output_polfile
READ（12,*）output_polfile
READ（12,*）output_polfile
READ（12,*）output_polfile
READ（12,*）output_polfile
READ（12,*）plane
READ（12,*）output_polfile
READ（12,*）sourcefile
READ（12,*）output_polfile
READ（12,*）stationfile
READ（12,*）output_polfile
READ（12,*）output_grafile
READ（12,*）output_polfile
READ（12,*）output_magfile
READ（12,*）output_polfile
CLOSE（12）
ENDSUBROUTINE
!
检查输入参数是否正确
SUBROUTINECHECK_CMD（component,order,plane,stationfile）
IMPLICITNONE
EXTERNALUPPER
INTEGERcomponent,order,plane,cc
CHARACTER*80stationfile,UPPER
CHARACTERsign
PRINT*,'请确认以下信息是否正确：
'
SELECTCASE（component）
CASE
（1）
PRINT*,'1计算重力或磁力异常'
CASE
（2）
PRINT*,'1计算重力或磁力异常的导数'
CASEDEFAULT
PRINT*,'错误：
场分量类型不合法，请修改参数文件！
！
'
STOP
ENDSELECT
IF（component==2.AND.order<0）THEN
PRINT*,'错误：
导数阶次不合法，请修改参数文件！
！
'
STOP
ENDIF
SELECTCASE（plane）
CASE
（1）
PRINT*,'2平面规则网'
CASE
（2）
PRINT*,'2平面非规则网'
CASE（3）
PRINT*,'2曲面规则网'
CASE（4）
PRINT*,'2曲面非规则网'
CASEDEFAULT
PRINT*,'错误：
测网类型不合法，请修改参数文件！
！
'
STOP
ENDSELECT
cc=LEN_TRIM（stationfile）
IF（（plane==1.OR.plane==3）.AND.
&（UPPER（stationfile（cc-3:
cc））.ne.'.GRD'））THEN
PRINT*,'错误：
测点文件类型不合法，请修改参数文件！
！
'
STOP
ENDIF
IF（（plane==2.OR.plane==4）.AND.
&（UPPER（stationfile（cc-3:
cc））.ne.'.TXT'））THEN
PRINT*,'错误：
测点文件类型不合法，请修改参数文件！
！
'
STOP
ENDIF
PRINT*,'以上内容是否正确？
（Y/N）'
READ*,sign
IF（UPPER（sign）=='N'）THEN
PRINT*,'请重新修改参数文件！
'
STOP
ENDIF
ENDSUBROUTINE
!
字符的小写转大写
FUNCTIONUPPER（string）
IMPLICITNONE
CHARACTER*（*）string,UPPER
INTEGERl,i,ASC
l=LEN（string）
DOi=1,l
ASC=ICHAR（string（i:
i））
IF（ASC>=97.AND.ASC<=122）THEN
ASC=ASC-32
string（i:
i）=CHAR（ASC）
ENDIF
ENDDO
UPPER=string
RETURN
ENDFUNCTION
!
读取场源的个数
SUBROUTINEGET_NUMBER_SOURCE（filename,n）
IMPLICITNONE
INTEGERn
CHARACTER*80filename
CHARACTER*20aa
OPEN（101,file=filename,action='read'）
n=0
DO
READ（101,*,END=100,ERR=100）aa
n=n+1
ENDDO
100CLOSE（101）
ENDSUBROUTINE
!
读取场源参数
SUBROUTINEINPUT_SOURCE（sourcefile,SOURCE,m,n）
IMPLICITNONE
INTEGERm,n,i,j
REALax,ay,incline
CHARACTER*80sourcefile
REALSOURCE（1:
m,1:
n）
OPEN（12,FILE=sourcefile,ACTION='READ'）
READ（12,*）（（SOURCE（i,j）,j=1,n）,i=1,m）
CLOSE（12）
DOi=1,m
ax=SOURCE（i,4）;ay=SOURCE（i,5）;incline=SOURCE（i,3）
IF（（ax+ay）==90..OR.（ax+ay）==270..OR.ABS（ax-ay）==90.）THEN
ELSE
PRINT*,'错误：
磁化方向参数不正确！
请修改场源参数文件！
！
'
PRINT*,'出错行：
',i
STOP
ENDIF
IF（incline<-90.OR.incline>90）THEN
PRINT*,'错误：
磁倾角参数不正确！
请修改场源参数文件！
！
'
PRINT*,'出错行：
',i
STOP
ENDIF
ENDDO
ENDSUBROUTINE
!
读取测量点的个数
SUBROUTINEGET_n_station（stationfile,plane,xmin,xmax,ymin,ymax,
&mpoint,line,n_station）
IMPLICITNONE
INTEGERplane,mpoint,line,n_station
REALxmin,xmax,ymin,ymax
CHARACTER*80stationfile
IF（plane==1.OR.plane==3）THEN
CALLGET_GRID（stationfile,mpoint,line,xmin,xmax,ymin,ymax）
n_station=mpoint*line
ELSE
mpoint=0;line=0;xmin=0;xmax=0;ymin=0;ymax=0
CALLGET_NUMBER_SOURCE（stationfile,n_station）
ENDIF
ENDSUBROUTINE
!
读取规则网测网网格参数
SUBROUTINEGET_GRID（inputfile,mpoint,line,xmin,xmax,ymin,ymax）
IMPLICITNONE
CHARACTER*80inputfile
INTEGERmpoint,line
REALxmin,xmax,ymin,ymax
OPEN（101,FILE=inputfile,ACTION='READ'）
READ（101,*）
READ（101,*）mpoint,line,xmin,xmax,ymin,ymax
CLOSE（101）
ENDSUBROUTINE
!
读取测点坐标
SUBROUTINEINPUT_STATION（stationfile,plane,xmin,xmax,ymin,ymax,
&mpoint,line,XYZ,n_station）
IMPLICITNONE
REALxmin,xmax,ymin,ymax,dx,dy,temp
INTEGERplane,mpoint,line,n_station
REALXYZ（1:
n_station,1:
3）
CHARACTER*80stationfile
REALi,j
OPEN（12,FILE=stationfile,ACTION='READ'）
IF（plane==1.OR.plane==3）THEN
DOi=1,5,1
READ（12,*）
ENDDO
dx=（xmax-xmin）/（mpoint-1）
dy=（ymax-ymin）/（line-1）
DOj=1,line
DOi=1,mpoint
XYZ（（j-1）*mpoint+i,1）=（i-1）*dx+xmin
XYZ（（j-1）*mpoint+i,2）=（j-1）*dy+ymin
ENDDO
ENDDO
IF（plane==1）THEN
READ（12,*）XYZ（1,3）
XYZ（:
3）=XYZ（1,3）
ELSE
READ（12,*）XYZ（:
3）
ENDIF
ELSEIF（plane==2.OR.plane==4）THEN
IF（plane==2）THEN
READ（12,*）XYZ（1,:
）
XYZ（:
3）=XYZ（1,3）
DOi=2,n_station
READ（12,*）XYZ（i,1）,XYZ（i,2）
ENDDO
ELSE
DOi=1,n_station
READ（12,*）XYZ（i,:
）
ENDDO
ENDIF
ELSE
PRINT*,'测网类型不合法！
'
ENDIF
ENDSUBROUTINE
!
计算各测点的异常值
SUBROUTINECAL_ANOM（SOURCE,n_source,cs,XYZ,n_station,cxyz,
&ANOM_GRA,ANOM_MAG,ANOM_POL,component）
IMPLICITNONE
EXTERNALDg,ElementDg,DT
REALDg,ElementDg,dg0,DT,dt0
INTEGERn_source,cs,n_station,cxyz,i,j,component
REALSOURCE（1:
n_source,1:
cs）,XYZ（1:
n_station,1:
cxyz）
REALANOM_GRA（1:
n_station）,ANOM_MAG（1:
n_station）
REALANOM_POL（1:
n_station）
ANOM_GRA=0.;ANOM_MAG=0.;ANOM_POL=0.
IF（component==1）THEN
DOi=1,n_station,1
DOj=1,n_source
dg0=Dg（XYZ（i,1）*1E6,XYZ（i,2）*1E6,XYZ（i,3）*1E6,SOURCE（j,6:
7）*1E6
&,SOURCE（j,8:
9）*1E6,SOURCE（j,10:
11）*1E6,SOURCE（j,1）,
&ElementDg）
ANOM_GRA（i）=ANOM_GRA（i）+dg0
dt0=DT（XYZ（i,1）*1E6,XYZ（i,2）*1E6,XYZ（i,3）*1E6,SOURCE（j,6:
7）*1E6
&,SOURCE（j,8:
9）*1E6,SOURCE（j,10:
11）*1E6,SOURCE（j,2）*1E-6,
&SOURCE（j,3）,SOURCE（j,4）,SOURCE（j,5））
ANOM_MAG（i）=ANOM_MAG（i）+dt0
dt0=DT（XYZ（i,1）*1E6,XYZ（i,2）*1E6,XYZ（i,3）*1E6,SOURCE（j,6:
7）*1E6
&,SOURCE（j,8:
9）*1E6,SOURCE（j,10:
11）*1E6,SOURCE（j,2）*1E-6,
&90.,SOURCE（j,4）,SOURCE（j,5））
ANOM_POL（i）=ANOM_POL（i）+dt0
ENDDO
ENDDO
ELSEIF（component==2）THEN
PRINT*,'重磁异常导数计算程序正在开发中...'
PRINT*,'$目前暂不支持此项功能$'
STOP
ENDIF
ENDSUBROUTINE
!
*******计算某源点在计算点引起的重力异常程序集开始*************
FUNCTIONDg（X,Y,Z,X_Source,Y_Source,Z_Source,Density,ElementDg）
IMPLICITNONE
REALX,Y,Z,X_Source（1:
2）,Y_Source（1:
2）,Z_Source（1:
2）,
&Density,ElementDg
REALDg,ELMT
INTEGERI,J,K,SIGN
Dg=0
ELMT=0
DOI=1,2
DOJ=1,2
DOK=1,2
SIGN=（-1）**（I+J+K）
ELMT=ElementDg（X,Y,Z,X_Source（I）,
&Y_Source（J）,Z_Source（K）,Density）
Dg=Dg+ELMT*SIGN
ENDDO
ENDDO
ENDDO
ENDFUNCTION
!
计算积分表达式中的某一分量
FUNCTIONElementDg（X,Y,Z,X_Source,Y_Source,Z_Source,Density）
IMPLICITNONE
REALG
PARAMETER（G=6.672E-8）
REALElementDg,X,Y,Z,X_Source,Y_Source,Z_Source,Density,R,T1,T2,T3
R=SQRT（（X_Source-X）**2+（Y_Source-Y）**2+（Z_Source-Z）**2）
T1=（X_Source-X）*LOG（（Y_Source-Y）+R）
T2=（Y_Source-Y）*LOG（（X_Source-X）+R）
T3=（Z_Source-Z）*ATAN（（X_Source-X）*（Y_Source-Y）/（（Z_Source-Z）*R））
ElementDg=-G*Density*（T1+T2-T3）
ENDFUNCTION
!
*******计算某源点在计算点引起的重力异常程序集结束*************
!
*******计算某源点在计算点引起的重力异常程序集开始*************
FUNCTIONDT（x,y,z,X_Source,Y_Source,Z_Source,moment,incline,ax,ay）
IMPLICITNONE
REALx,y,z,moment,incline,ax,ay,DT
REALX_Source（1:
2）,Y_Source（1:
2）,Z_Source（1:
2）
INTEGERi,j,k,sign
REALhax,hay,za,mx,my,mz,pi,r,xs,ys,zs
pi=3.14159265358
mx=moment*COSD（incline）*COSD（ax）
my=moment*COSD（incline）*COSD（ay）
mz=moment*SIND（incline）
hax=0.;hay=0.;za=0.
DOi=1,2,1
DOj=1,2,1
DOk=1,2,1
sign=（-1）**（i+j+k）
xs=X_Source（i）;ys=Y_Source（j）;zs=Z_Source（k）
r=SQRT（（x-xs）**2+（y-ys）**2+（z-zs）**2）
hax=hax+sign*（-mx*ATAN（（xs-x）*（ys-y）/（（xs-x）**2+r*（zs-z）+
&（zs-z）**2））+my*LOG（r+zs-z）+mz*LOG（r+ys-y））
hay=hay+sign*（mx*LOG（r+zs-z）-my*ATAN（（xs-x）*（ys-y）/（
&（ys-y）**2+r*（zs-z）+（zs-z）**2））+mz*LOG（r+xs-x））
za=za+sign*（mx*LOG（r+ys-y）+my*LOG（r+xs-x）-mz*
&ATAN（（xs-x）*（ys-y）/r/（zs-z）））
ENDDO
ENDDO
ENDDO
DT=hax*COSD（incline）*COSD（ax）+hay*COSD（incline）*COSD（ay）+
&za*SIND（incline）
DT=DT/4/pi
ENDFUNCTION
!
*******计算某源点在计算点引起的重力异常程序集结束*************
!
输出重磁异常到外部文件中
SUBROUTINEOUTPUT（plane,XYZ,ANOM_GRA,ANOM_MAG,ANOM_POL,n_station
&,xmin,xmax,ymin,ymax,mpoint,line,output_grafile,output_magfile
&,output_polfile）
INTEGERplane,n_station,mpoint,line,i,j
REALANOM_GRA（1:
n_station）,ANOM_MAG（1:
n_station）
REALANOM_POL（1:
n_station）,XYZ（1:
n_station,1:
3）
REALxmin,xmax,ymin,ymax,gzmin,gzmax,mzmin,mzmax,pzmin,pzmax
REALeigv
CHARACTER*80output_grafile,output_magfile,output_polfile
eigv=1.7014100091878E+038
OPEN（12,FILE=output_grafile,ACTION='WRITE'）
OPEN（13,FILE=output_magfile,ACTION='WRITE'）
OPEN（14,FILE=output_polfile,ACTION='WRITE'）
IF（plane==1.OR.plane==3）THEN
gzmin=eigv;gzmax=-eigv;mzmin=eigv
mzmax=-eigv;pzmin=eigv;pzmax=-eigv
DOi=1,n_station,1
IF（ANOM_G
升级会员 