糖尿病康复 > GRACE同震重力

GRACE同震重力

时间：2023-10-25 05:49:19

相关推荐

GRACE同震重力

!*****************************************************************************************************

!—————————————————————————————————————

! 本次程序旨在写出一套计算同震形变的程序，为将来的高分辨率的卫星重力用于地震研究打下基础

!—————————————————————————————————————

! 功能：时变重力场计算同震形变

!—————————————————————————————————————

! 程序说明：此程序主要包括时变重力场的重力计算，重力扰动计算，同震重力扰动的提取，断层参数的非线性反演

! 断层滑动的线性反演，以及对提取的同震形变的验证

!—————————————————————————————————————

! 参数介绍：

! n:重力场阶数（此次直接采用60阶的时变重力场）

! m:重力场次数

! G:万有引力常数

! Snm:球谐系数

! Pnm:规格化的勒让德函数

! R:地球平均半径

! r:计算点到地球质心的距离

! M1:地球质量

! gravity(:,:):计算点的重力扰动值

! 输入文件：GSM-2_032-059_GRAC_UTCSR_BA01_0600，一定时间段内的时变重力场系数

! 输出文件：output_file_field.grd，计算得到的一定时间段的重力扰动

!—————————————————————————————————————-

!******************************************************************************************************

program shibian_codeformation

implicit none

real:: R=6.371*1E6,G=6.67*1E-11 !M,N·m²/kg²

real*8 Cnm(0:60,0:60),Snm(0:60,0:60)

real*8 Pnm(0:61,0:61,0:90)!因为余纬的0到90度，0.5度一个单位

real*8 d_g(0:90,0:180)!时变重力

real*8 Cnm1(0:60,0:60),Snm1(0:60,0:60)

REAL:: M1=5.965*1E24 !kg

!real*8 average_field(1:64)

real*4 Rnm(0:61,0:61,0:90,0:180),SSnm(0:61,0:61,0:90,0:180)!60阶的时变重力场的分辨率为3度

integer i,j

!计算完全规格化缔和勒让德函数

call read_parameter1(Pnm)

!PAUSE

!读取球谐系数

call input_Cnm_parameter(Cnm,Snm,Cnm1,Snm1)

!读取代换的C20

call readparameterc20(Cnm)

Cnm=Cnm-Cnm1

Snm=Snm-Snm1

!计算时变重力场的变化

call cogravity(G,M1,R,Cnm,Snm,Pnm,d_g,Rnm,SSnm)

!计算高斯和去各向同性之后的重力场变化

!call cogravity_Gauss_flter(G,M1,R,Cnm,Snm,Pnm,d_g,Rnm,SSnm)

!输出计算所得到的重力场的变化

call output_cogravity(d_g)

!输出与平均重力值差异的重力场变化

!call average_output_cogravity(average_field)

end program

!******输入部分**********************************************************

!—————————————————-

!功能:从参数文件中读取有关参数,利用EGM作为背景场，计算得到了时变重力场与EGM之间系数的差别

! 输入:时变重力场文件,参看文献24,

! 输出:Cnm,Snm

!—————————————————-

subroutine input_Cnm_parameter(Cnm,Snm,Cnm1,Snm1)

implicit none

real*8 Cnm(0:60,0:60),Snm(0:60,0:60)

real*8 Cnm1(0:60,0:60),Snm1(0:60,0:60)

integer i,j,k,kk,l

integer n,m

character a,f

real b,c,d,e

real*8 jiecheng1,jiecheng2,dd

Cnm=0.0

Snm=0.0

Cnm1=0.0

Snm1=0.0

open(10,file=’GSM-2_032-059_GRAC_UTCSR_BA01_0600’)

do k=1,1891,1

read(10,*)a,n,m,Cnm(n,m),Snm(n,m),b,c,d,e,f

enddo

close(10)

open(10,file=’GSM-2_032-059_GRAC_UTCSR_BA01_0600’)

do k=1,1891,1

read(10,*)a,n,m,Cnm1(n,m),Snm1(n,m),b,c,d,e,f

enddo

close(10)

!!计算完全规格化的球谐系数

!do i=0,60,1

! do j=0,i,1

! k=i-j

! kk=i+j

! if(k==0)then

! jiecheng1=1.0

! else

! do l=1,k,1

! jiecheng1=jiecheng1*l

! enddo

! endif

! if(kk==0)then

! jiecheng2=1.0

! else

! do l=1,kk,1

! jiecheng2=jiecheng2*kk

! enddo

! endif

!! !公式中的K的计算

! if(j==0)then

! dd=1

! else

! dd=0

! endif

! Cnm(i,j)=sqrt(ABS((dd*(2.0*i+1)*jiecheng1)/jiecheng2))*Cnm(i,j)

! Snm(i,j)=sqrt(ABS((dd*(2.0*i+1)*jiecheng1)/jiecheng2))*Snm(i,j)

!!!

!! Cnm1(i,j)=sqrt(ABS((dd*(2.0*i+1)*jiecheng1)/jiecheng2))*Cnm1(i,j)

!! Snm1(i,j)=sqrt(ABS((dd*(2.0*i+1)*jiecheng1)/jiecheng2))*Snm1(i,j)

!!!

! jiecheng1=1.0

! jiecheng2=1.0

! enddo

!enddo

end subroutine

!—————————————————-

!功能:从参数文件中读取C20有关参数

! 输入:时变重力场文件

! 输出:C20

!—————————————————-

subroutine readparameterc20(Cnm)

implicit none

real*8 Cnm(0:60,0:60)

real a,b,d,e

integer n,m

open(10,file=’SLR_C20.txt’,status=’old’)

read(10,*)

!read(10,*)

read(10,*)a,b,Cnm(2,0),d,e

close(10)

!Cnm(2,0)=SQRT(ABS(2.0/10.0))!替换C20

end subroutine

!—————————————————-

!功能:计算规格化勒让德函数：参看文献23

!输入:

!输出:Pnm

!—————————————————-

subroutine read_parameter1(P)

real*8 P(0:61,0:61,0:90)

integer i,j,k,kk,l

real*8 k1

real*8 anm,bnm

real*8 Nnm(0:60,0:60)!完全规格化系数

real*8:: jiecheng1=1.0,jiecheng2=1.0

real dd

KK=1

Nnm=0.0

Pnm=0.0

anm=0.0

bnm=0.0

do i=0,60,1

do j=0,i,1

do k=0,90,1

k1=k*2.0/90*3.14!3度乘回去，将经度角度制转化为弧度制!参考文献23算法2if(j==0)thenp(0,0,k)=1.0p(1,0,k)=sqrt(3.0)*cos(k1) if(i>1)thenp(i,0,k)=sqrt(4.0*i*i-1.0)/i*cos(k1)*p(i-1,0,k)-(i-1)/i*&&sqrt((2.0*i+1.0)/(2*i-3.0))*p(i-2,0,k)!write(*,*)P(i,0,k)endifendif p(1,1,k)=sqrt(3.0)*sin(k1) if(j==i.and.(i>1.and.j>1))thenp(i,j,k)=sqrt((2.0*i+1)/(2*i))*sin(k1)*p(i-1,j-1,k)endifif(j==i-1)thenp(i,j,k)=sqrt(2.0*i+1)*sin(k1)*p(i-1,i-1,k)endifif(j<(i-1))thenp(i,j,k)=sqrt((4.0*i*i-1)/(i*i-j*j))*sin(k1)*p(i-1,j,k)&&-sqrt(((2*i+1)*((i-1)**2-j**2))/((2*i-3.0)*(i**2-j**2)))*P(i-2,j,k)endif !write(*,*)p(i,j,k),i,jenddoenddo

enddo

!k=0

!kk=0

!!计算完全规格化系数

!do i=0,60,1

! do j=0,i,1

! if(j==0)then

! dd=1

! else

! dd=0

! endif

! k=i-j

! kk=i+j

! if(k==0)then

! jiecheng1=1.0

! else

! do l=1,k,1

! jiecheng1=jiecheng1*l

! enddo

! endif

! if(kk==0)then

! jiecheng2=1.0

! else

! do l=1,kk,1

! jiecheng2=jiecheng2*kk

! enddo

! endif

! Nnm(i,j)=sqrt(ABS((2.0-dd)*(2*i+1)*jiecheng1/jiecheng2))

! !Nnm(i,j)=sqrt(ABS(1.0))

! !write(,)Nnm(0,0),k,kk

! !PAUSE

! jiecheng1=1.0

! jiecheng2=1.0

! enddo

!enddo

!!pause

!!计算完全规格化缔和勒让德函数

!do i=0,60,1

! do j=0,i,1

! do k=0,90,1

! P(i,j,k)=Nnm(i,j)*P(i,j,k)

! !Write(,)P(i,j,k),i,j,k

! enddo

!enddo

!pause

!open (10,file=’X.txt’)

!do i=0,20,1

! do j=0,i,1

! do k=0,90,1

!! ! do l=561,601,1

! write(10,*)i,j,k,P(i,j,k)

! ! enddo

! enddo

!enddo

!CLOSE(10)

!pause

end subroutine

!—————————————————-

!功能:计算时变重力扰动（此时认为R=r),参考文献11

!输入:Snm,Cnm,Pnm

!输出:cogrvity

!—————————————————-

subroutine cogravity(G,M,R,Cnm,Snm,Pnm,d_g,Rnm,SSnm)

implicit none

real G,M,R,C

real*8 Cnm(0:60,0:60),Snm(0:60,0:60)

real*8 Pnm(0:61,0:61,0:90)

real*4 Rnm1(0:61,0:61,0:90,0:180),SSnm1(0:61,0:61,0:90,0:180)&

&,Rnm(0:61,0:61,0:90,0:180),SSnm(0:61,0:61,0:90,0:180)!60阶的时变重力场的分辨率为3度

real*8 d_g(0:90,0:180)

real*8 f1,f2

integer i,j,k,l

real k1

real*8 jiecheng

C=G*M/(R**2)*1E8 !GM/R**2常数值,微伽

!write(,)C

!pause

Rnm=0.0

SSnm=0.0

Rnm1=0.0

SSnm1=0.0

!计算Rnm,SSnm的具体值

do i=0,60,1

do j=0,i,1

do k=0,90,1

do l=0,180,1

k1=(j*1*2)/90*3.141592653

Rnm1(i,j,k,l)=Pnm(i,j,k)*cos(k1)

SSnm1(i,j,k,l)=Pnm(i,j,k)*sin(k1)

Rnm(i,j,k,l)=jiecheng(i,j)*Rnm1(i,j,k,l)

SSnm(i,j,k,l)=jiecheng(i,j)*SSnm1(i,j,k,l)

!write(,) Pnm(i,j,k), SSnm(i,j,k,l)

enddo

!pause

enddo

!write(,) jiecheng(i,j)

enddo

!pause

d_g=0.0

f1=0.0

f2=0.0

!计算时变重力值,先余纬再经度

do i=0,90,1

do j=0,180,1

do k=2,60,1

do l=0,k,1

f1=f1+Cnm(k,l)*Rnm(k,l,i,j)+Snm(k,l)*Ssnm(k,l,i,j)

enddof2=f2+(k+1)*f1f1=0.0!PAUSEenddod_g(i,j)=C*f2!write(*,*) f2,d_g(i,j)f2=0.0

enddo

!pause

end subroutine

!—————————————————-

!功能:计算时变重力变化值公式中的阶乘

!输入:

!输出:

!—————————————————-

real*8 function jiecheng(n,m)

implicit none

integer n,m

real*8 ji1,ji2

integer i,j,k,l

i=n-m

j=n+m

ji1=1.0

ji2=1.0

if(i==0)then

ji1=1.0

else

do k=1,i,1

ji1=ji1*k

enddo

endif

if(j==0)then

ji2=1.0

else

do l=1,j,1

ji2=ji2*l

enddo

endif

!write(,)ji1/ji2,n

if (m==0)then

jiecheng=sqrt((2.0*n+1)*(ji1/ji2))

else

jiecheng=sqrt(2.0*(2*n+1)*(ji1/ji2))

endif

end function

!—————————————————-

!功能:对计算得到的卫星重力异常进行滤波（各向同性高斯滤波，去相关滤波)选择滤波半径是300KM

!参考文献22

!输入:

!输出:

!—————————————————-

subroutine cogravity_Gauss_flter(G,M,R,Cnm,Snm,Pnm,d_g,Rnm,SSnm)

implicit none

real*8 d_g(0:90,0:180)

real G,M,R,C

real*8 Cnm(0:60,0:60),Snm(0:60,0:60)

real*8 Pnm(0:61,0:61,0:61)

real*4 Rnm(0:61,0:61,0:90,0:180),SSnm(0:61,0:61,0:90,0:180)

real*8 f1,f2

integer i,j,k,l

real*8 jiecheng

real*8 W(0:60)

real*8:: r1=300,alpha!滤波因子，滤波半径

C=G*M/(R**2)*1E8 !GM/R**2常数值,微伽

!计算滤波因子

W(0)=1.0

alpha=ALOG(2.0)/(1.0-cos(r1/6371))

W(1)=(1+EXP(-2.0*alpha))/(1-EXP(-2.0*alpha))-1.0/alpha

do i=2,60,1

W(i)=-(2.0*i-1.0)/(alpha)*W(i-1)+W(i-2)

enddo

f1=0.0

f2=0.0

!计算各向同性高斯滤波之后的时变重力值,先余纬再经度

do i=0,90,1

do j=0,180,1

do k=2,60,1

do l=0,k,1

f1=f1+Cnm(k,l)*Rnm(k,l,i,j)+Snm(k,l)*Ssnm(k,l,i,j)

enddo

f2=f2+(k+1)*f1*W(k)

f1=0.0

enddo

d_g(i,j)=C*f2

!write(,) f2,d_g(i,j)

f2=0.0

enddo

!PAUSE

!进行去相关化处理

end subroutine

!*********输出部分**********************************************************

!—————————————————-

!功能:输出计算所得到的重力变化

!输入:d_g

!输出:d_g.grad文件

!—————————————————-

subroutine output_cogravity(d_g)

implicit none

real min,max

real*8 d_g(0:90,0:180)

real:: xmin,xmax,ymin,ymax

integer i,j,k,n,m,xy

real*8 field(1:16471)!计算点总共有121个点

real zmin,zmax

m=91

n=181

!write(,)d_g

open(20,file=’output.2.grd’)

ymax=90

ymin=-90

xmin=0

xmax=360

do j=0,90,1

do k=0,180,1

field(j*91+k+1)=d_g(j,k)

!WRITE(,) d_g(j,k)

enddo

!open (10,file=’X.txt’)

!do i=0,90,1

! do j=0,180,1

! write(10,*)i,j,d_g(i,j)

! enddo

!enddo

!CLOSE(10)

!!pause

!PAUSE

zmin=field(1)

zmax=field(1)

do i= 2,16471,1

zmin=min(zmin,field(i))

zmax=max(zmax,field(i))

end do

write (20,’(A)’)’DSAA’

write (20,*)n,m

write (20,*)xmin,xmax

write (20,*)ymin,ymax

write (20,*)zmin,zmax

do i=0,180,1

xy=i*91

write(20,*) (field(xy+j),j=1,91,1)

enddo

close(20)

end subroutine

!—————————————————-

!功能:输出计算所得到的平均重力变化

!输入:

!输出:所有月份的平均重力变化文件和所有月份与的差值变化文件

!—————————————————-

subroutine average_output_cogravity(average_field)

implicit none

real min,max

real*8 field2(1:64),field3(1:64),a(1:64)

!real*8 field4(1:64),field5(1:64),field6(1:64),&

! &field7(1:64),field8(1:64),field9(1:64),field10(1:64),field11(1:64),field12(1:64),field01(1:64),field02(1:64),&

! &field03(1:64),field04(1:64),field05(1:64),field06(1:64),field07(1:64),field08(1:64),field09(1:64),field010(1:64),&

! &field011(1:64),field012(1:64),field112(1:64)

real::xmin,xmax,ymin,ymax

integer i,j,k,n,m,xy

real*8 average_field(1:64)!计算点总共有121个点

real zmin,zmax

!real*8 field24(1:24,1:64)

!读取所有月份的时变重力值

open(20,file=’output.2.grd’,status=’old’)

read (20,*)

read (20,*)m,n

read (20,*)xmin,xmax

read (20,*)ymin,ymax

read (20,*)zmin,zmax

do i=33,40,1

xy=(i-33)*8

read(20,*) (field2(xy+j),j=1,8,1)

enddo

close(20)

open(20,file=’output.31.grd’,status=’old’)

read (20,*)

read (20,*)m,n

read (20,*)xmin,xmax

read (20,*)ymin,ymax

read (20,*)zmin,zmax

do i=33,40,1

xy=(i-33)*8

read(20,*) (field3(xy+j),j=1,8,1)

enddo

close(20)

!open(20,file=’output.4.grd’,status=’old’)

!read (20,*)

!read (20,*)m,n

!read (20,*)xmin,xmax

!read (20,*)ymin,ymax

!read (20,*)zmin,zmax

!do i=33,40,1

! xy=(i-33)*8

! read(20,*) (field4(xy+j),j=1,8,1)

!enddo

!close(20)

!open(20,file=’output.5.grd’,status=’old’)

!read (20,*)

!read (20,*)m,n

!read (20,*)xmin,xmax

!read (20,*)ymin,ymax

!read (20,*)zmin,zmax

!do i=33,40,1

! xy=(i-33)*8

! read(20,*) (field5(xy+j),j=1,8,1)

!enddo

!close(20)

!open(20,file=’output.6.grd’,status=’old’)

!read (20,*)

!do i=33,40,1

! xy=(i-33)*8

! read(20,*) (field6(xy+j),j=1,8,1)

!enddo

!close(20)

!open(20,file=’output.7.grd’,status=’old’)

!read (20,*)

!do i=33,40,1

! xy=(i-33)*8

! read(20,*) (field7(xy+j),j=1,8,1)

!enddo

!close(20)

!open(20,file=’output.8.grd’,status=’old’)

!read (20,*)

!do i=33,40,1

! xy=(i-33)*8

! read(20,*) (field8(xy+j),j=1,8,1)

!enddo

!close(20)

!open(20,file=’output.9.grd’,status=’old’)

!read (20,*)

!do i=33,40,1

! xy=(i-33)*8

! read(20,*) (field9(xy+j),j=1,8,1)

!enddo

!close(20)

!open(20,file=’output.10.grd’,status=’old’)

!read (20,*)

!do i=33,40,1

! xy=(i-33)*8

! read(20,*) (field10(xy+j),j=1,8,1)

!enddo

!close(20)

!open(20,file=’output.11.grd’,status=’old’)

!read (20,*)

!do i=33,40,1

! xy=(i-33)*8

! read(20,*) (field11(xy+j),j=1,8,1)

!enddo

!close(20)

!open(20,file=’output.12.grd’,status=’old’)

!read (20,*)

!do i=33,40,1

! xy=(i-33)*8

! read(20,*) (field12(xy+j),j=1,8,1)

!enddo

!close(20)

!open(20,file=’output.1.grd’,status=’old’)

!read (20,*)

!do i=33,40,1

! xy=(i-33)*8

! read(20,*) (field01(xy+j),j=1,8,1)

!enddo

!close(20)

!open(20,file=’output.2.grd’,status=’old’)

!read (20,*)

!do i=33,40,1

! xy=(i-33)*8

! read(20,*) (field02(xy+j),j=1,8,1)

!enddo

!close(20)

!open(20,file=’output.3.grd’,status=’old’)

!read (20,*)

!do i=33,40,1

! xy=(i-33)*8

! read(20,*) (field03(xy+j),j=1,8,1)

!enddo

!close(20)

!open(20,file=’output.4.grd’,status=’old’)

!read (20,*)

!do i=33,40,1

! xy=(i-33)*8

! read(20,*) (field04(xy+j),j=1,8,1)

!enddo

!close(20)

!open(20,file=’output.5.grd’,status=’old’)

!read (20,*)

!do i=33,40,1

! xy=(i-33)*8

! read(20,*) (field05(xy+j),j=1,8,1)

!enddo

!close(20)

!open(20,file=’output.6.grd’,status=’old’)

!read (20,*)

!do i=33,40,1

! xy=(i-33)*8

! read(20,*) (field06(xy+j),j=1,8,1)

!enddo

!close(20)

!open(20,file=’output.7.grd’,status=’old’)

!read (20,*)

!do i=33,40,1

! xy=(i-33)*8

! read(20,*) (field07(xy+j),j=1,8,1)

!enddo

!close(20)

!open(20,file=’output.8.grd’,status=’old’)

!read (20,*)

!do i=33,40,1

! xy=(i-33)*8

! read(20,*) (field08(xy+j),j=1,8,1)

!enddo

!close(20)

!open(20,file=’output.9.grd’,status=’old’)

!read (20,*)

!do i=33,40,1

! xy=(i-33)*8

! read(20,*) (field09(xy+j),j=1,8,1)

!enddo

!close(20)

!open(20,file=’output.10.grd’,status=’old’)

!read (20,*)

!do i=33,40,1

! xy=(i-33)*8

! read(20,*) (field010(xy+j),j=1,8,1)

!enddo

!close(20)

!open(20,file=’output.11.grd’,status=’old’)

!read (20,*)

!do i=33,40,1

! xy=(i-33)*8

! read(20,*) (field011(xy+j),j=1,8,1)

!enddo

!close(20)

!open(20,file=’output.12.grd’,status=’old’)

!read (20,*)

!do i=33,40,1

! xy=(i-33)*8

! read(20,*) (field012(xy+j),j=1,8,1)

!enddo

!close(20)

!open(20,file=’output.2.grd’,status=’old’)

!read (20,*)

!do i=33,40,1

! xy=(i-33)*8

! read(20,*) (field112(xy+j),j=1,8,1)

!enddo

!close(20)

!!!求平均值

!do i=64,1

! average_field(i)=1.0/24*(field2(i)+field3(i)+field4(i)+field5(i)+field6(i)+field7(i)+field8(i)+field9(i)+field10(i)+&

! &field11(i)+field12(i)+field01(i)+field02(i)+field03(i)+field04(i)+field05(i)+field06(i)+field07(i)+field08(i)+field09(i)+&

! &field010(i)+field011(i)+field012(i)+field112(i))

!enddo

!!输出平均值

!open(20,file=’outputaverage.grd’,status=’unknown’)

!write (20,’(A)’)’DSAA’

!write (20,*)m,n

!write (20,*)xmin,xmax

!write (20,*)ymin,ymax

!write (20,*)zmin,zmax

!do i=33,40,1

! xy=(i-33)*8

! write(20,*) (average_field(xy+j),j=1,8,1)

!enddo

!close(20)

!do j=1,24,1

! field24(1,j)=field2(j)

! field24(2,j)=field3(j)

! field24(3,j)=field4(j)

! field24(4,j)=field5(j)

! field24(5,j)=field6(j)

! field24(6,j)=field7(j)

! field24(7,j)=field8(j)

! field24(8,j)=field9(j)

! field24(9,j)=field10(j)

! field24(10,j)=field11(j)

! field24(11,j)=field12(j)

! field24(12,j)=field01(j)

! field24(13,j)=field02(j)

! field24(14,j)=field03(j)

! field24(15,j)=field04(j)

! field24(16,j)=field05(j)

! field24(17,j)=field06(j)

! field24(18,j)=field07(j)

! field24(19,j)=field08(j)

! field24(20,j)=field09(j)

! field24(21,j)=field010(j)

! field24(22,j)=field011(j)

! field24(23,j)=field012(j)

! field24(24,j)=field112(j)

!enddo

!open(20,file=’output_bianhua.grd’)

!write (20,’(A)’)’DSAA’

!write (20,*)m,n

!write (20,*)xmin,xmax

!write (20,*)ymin,ymax

!write (20,*)zmin,zmax

!!do k=1,24,1

!! do i=33,40,1

!! xy=(i-33)*8

!! write(20,*) (field24(k,(xy+j))-average_field(xy+j),j=1,8,1)

!! enddo

!! write(20,*)

!!enddo

!do i=1,64,1

! a(i)= field3(i)-field2(i)

!enddo

!do i=33,40,1

! xy=(i-33)*8

! write(20,*) (a(xy+j),j=1,8,1)

!enddo

!close(20)

end subroutine

如果觉得《GRACE同震重力》对你有帮助，请点赞、收藏，并留下你的观点哦！

本内容不代表本网观点和政治立场，如有侵犯你的权益请联系我们处理。

网友评论

网友评论仅供其表达个人看法，并不表明网站立场。