ROck Physics Toolbox  1.0 A microgeodynamics-based toolkit for rock physics.
main.f90
Go to the documentation of this file.
1 PROGRAM main
2  USE global
4  USE regional
5  !TYPE(unit_cell)::PREM
6  TYPE(region)::hawaii
7  TYPE(composition)::comp
8  TYPE(mantle)::model
9  INTEGER::ii
10
11  ! Define composition space
12  comp=composition(potential_temperature=1550,basalt_fraction=18)
14  ! Load the regional data
15  hawaii%NGRID=2
16
18  hawaii%MAGMA=set_melt(2) !
19  hawaii=melt_calculate(hawaii)
20
21  CALL vtk_write(hawaii,model)
22  print*,'Potential Temperature', comp%potential_temperature
23  print*,'basalt_fraction',comp%basalt_fraction
24  write(*,'(6(1x,a))'),'Depth to LVL',' Depth to 410',' Potential T',' Vs_obs', ' Vs_sol',' Melt Fraction'
25  DO ii=1,hawaii%NGRID
26  write(*,'(6(1x,f12.5))') hawaii%TOPO(ii,1),hawaii%TOPO(ii,2),hawaii%CELL(ii)%temperature,hawaii%CELL(ii)%vs_obs, hawaii%CELL(ii)%vs_sol, hawaii%CELL(ii)%melt_fraction
27  END DO
28  !CALL inversion_test(HAWAII)
29
30 CONTAINS
31
32  SUBROUTINE inversion_test(reg)
33  IMPLICIT NONE
34  TYPE(region),INTENT(in)::reg
35  TYPE(region)::reg2
36  INTEGER::ii
37  REAL(sp)::cont
38  REAL(sp),DIMENSION(4)::temp
39
40
41  reg2=reg
42
43  reg2%CELL(1:reg%NGRID)%melt_fraction=linspace_r(0.001_sp,0.1_sp,reg%NGRID)
44
45  reg2%CELL(1:reg%NGRID)%K_sol=655.6e9_sp!GPa
46  reg2%CELL(1:reg%NGRID)%G_sol=293.8e9_sp!GPa
47  reg2%CELL(1:reg%NGRID)%vp_sol=13716.6_sp! m/s
48  reg2%CELL(1:reg%NGRID)%vs_sol=7264.66_sp
49  reg2%CELL(1:reg%NGRID)%rho=5566.45_sp
50  reg2%CELL(1:reg%NGRID)%nu=0.3051_sp
51
52  reg2%CELL(1:reg%NGRID)%Kl=452.98_sp!GPa
53  reg2%CELL(1:reg%NGRID)%rhol=2.4_sp*2690.0_sp
54  print*,'================================================='
55  print*,'Error testing for melt fraction inversion'
56  print*,'================================================='
57  WRITE(*,'(5(1x,a))'),'|Original', ' |Nonlinear', '|Error', '| Bisection', '| Error|'
58  print*,'================================================='
59  DO ii=1,reg%NGRID
60  cont=whm(reg2%MAGMA%dihedral,reg2%CELL(ii)%melt_fraction)
61  temp=elastic_melt(cont,reg2%CELL(ii)%melt_fraction&
62  &,reg2%CELL(ii)%K_sol,reg2%CELL(ii)%G_sol,reg2%CELL(ii)%nu,&
63  &reg2%CELL(ii)%rho,reg2%CELL(ii)%Kl,reg2%CELL(ii)%rhol)
64  reg2%CELL(ii)%vs_eff=temp(1)
65  reg2%CELL(ii)%vp_eff=temp(2)
66  reg2%CELL(ii)%K_eff=temp(3)
67  reg2%CELL(ii)%G_eff=temp(4)
68  reg2%CELL(ii)%vs_obs=nonlinear_solve(0.0001_sp,0.45_sp,1.0e-4_sp,reg2%CELL(ii),reg2%MAGMA)
69  reg2%CELL(ii)%vp_obs=bisection(0.0001_sp,0.45_sp,1.0e-4_sp,reg2%CELL(ii),reg2%MAGMA)
70  WRITE(*,'(5(1x,es10.2))'),reg2%CELL(ii)%melt_fraction,reg2%CELL(ii)%vs_obs,abs(reg2%CELL(ii)%vs_obs-reg2%CELL(ii)%melt_fraction),reg2%CELL(ii)%vp_obs,abs(reg2%CELL(ii)%vp_obs-reg2%CELL(ii)%melt_fraction)
71  END DO
72  END SUBROUTINE inversion_test
73
74 END PROGRAM main