[Fortran] 纯文本查看 复制代码
PROGRAM GRAVITY
PARAMETER(N=50)
INTEGER XL
REAL X1,X2,X3,X4,X5 !层块状模型横坐标
REAL H11,H12,H13,H14,H15 !第一层顶点深度,KM
REAL H21,H22,H23,H24,H25 !第二层顶点深度,KM
REAL H31,H32,H33,H34,H35 !第三层顶点深度,KM
REAL GOBS0(N),GOBS1(N),GOBS2(N),GOBS3(N),GOBS4(N)
REAL GOBS5(N),GOBS6(N),GOBS7(N),GOBS8(N) !正演重力异常值
X1=2.1;X2=8;X3=14;X4=20;X5=22 !层块状模型横坐标
H11=1.0;H12=1.5;H13=0.5;H14=3.0;H15=1.0
H21=3.0;H22=2.5;H23=2.0;H24=4.0;H25=2.5
H31=4.5;H32=3.5;H33=4.0;H34=5.0;H35=4.0
RD1=0.1;RD2=0.2;RD3=0.1;RD4=0.4;RD5=0.4;RD6=0.1;RD7=0.2;RD8=0.1 !每个模块的剩余密度，单位g/cm3
DO I=0,50,1
CALL gobso(X1,X2,H11,H12,H21,H22,RD1) !把第一模块所有参数带进去
ENDDO
DO I=0,50,1
CALL gobso(X2,X3,H12,H13,H22,H23,RD2) !把第二模块所有参数带进去
ENDDO
DO I=0,50,1
CALL gobso(X3,X4,H13,H14,H23,H24,RD3) !把第三模块所有参数带进去
ENDDO
DO I=0,50,1
CALL gobso(X4,X4,H14,H15,H24,H25,RD4) !把第四模块所有参数带进去
ENDDO
DO I=0,50,1
CALL gobso(X1,X2,H21,H22,H31,H32,RD5) !把第五模块所有参数带进去
ENDDO
DO I=0,50,1
CALL gobso(X2,X3,H22,H23,H32,H33,RD6) !把第六模块所有参数带进去
ENDDO
DO I=0,50,1
CALL gobso(X3,X4,H23,H24,H33,H34,RD7) !把第七模块所有参数带进去
ENDDO
DO I=0,50,1
CALL gobso(X4,X5,H24,H25,H34,H35,RD8) !把第八模块所有参数带进去
ENDDO
DO I=0,50,1
GOBS0(I)=GOBS1(I)+GOBS2(I)+GOBS3(I)+GOBS4(I)+GOBS5(I)+GOBS6(I)+GOBS7(I)+GOBS8(I)
PRINT*,GOBS0(I) !
ENDDO
CONTAINS
SUBROUTINE gobso(X1,X2,H1,H2,H3,H4,RD)
INTEGER I
REAL H1,H2 !梯形上侧边顶点深度值,m
REAL H3,H4 !梯形下侧边顶点深度值,m
REAL X1,X2 !梯形上侧边横坐标
REAL X !测点位置
REAL RD !剩余密度
REAL A,B,C,L,P,K,B0A,BTA,gobs,AK
REAL AO,BO,CO,PO,KO,B0AO,BTAO,AKO
REAL EX,CX,EX0,CX0
REAL gobsOO,gobsOT,gobsTO,gobsTT
G=0.06672 !万有引力常数
L=1/(X2-X1)
K=L*(H2-H1)
KO=L*(H4-H3)
A=1+K**2
AO=1+KO**2
X=-11
DO I=1,50,1
X=X+1
P=K*(X-X1)+H1
PO=KO*(X-X1)+H3
B=2*K*P
B=2*KO*P1
C=P**2
C0=PO**2
EX=X1-X
CX=ABS(A*EX**2+B*EX+C)
EXO=X2-X
CXO=ABS(AO*EXO**2+BO*EXO+CO)
IF(P==0)THEN
BOA=0
BTA=0
ELSE
AK=ATAN((2*A*(X1-X)+B)/2*ABS(P))/(2*ABS(P))
B0A=(B**2-2*A*C)*AK/(2*A**2)
BTA=B*AK/(2*A**2)
ENDIF
IF(PO==O)THEN
BOAO=0
BTAO=0
ELSE
AKO=ATAN((2*AO*(X2-X)+BO)/2*ABS(PO))/(2*ABS(PO))
B0AO=(BO**2-2*AO*CO)*AKO/(2*AO**2)
BTAO=BO*AKO/(2*AO**2)
ENDIF
gobsOO=G*RD*(EX*LOG(CX)+2*A*(EX/A-B*LOG(CX)/(2*A)+BOA)-B*(LOG(CX)/(2*A)-BTA))
gobsOT=G*RD*(EX*LOG(CX)+2*AO*(EX/AO-BO*LOG(CX)/(2*AO)+BOAO)-BO*(LOG(CX)/(2*AO)-BTAO))
gobsTO=G*RD*(EXO*LOG(CXO)+2*A*(EXO/A-B*LOG(CXO)/(2*A)+BOA)-B*(LOG(CXO)/(2*A)-BTA))
gobsTT=G*RD*(EXO*LOG(CXO)+2*AO*(EXO/AO-BO*LOG(CXO)/(2*AO)+BOAO)-BO*(LOG(CXO)/(2*AO)-BTAO))
gobs=gobsTO+gobsOT-gobsOO-gobsTT
ENDDO
END SUBROUTINE
END