ROck Physics Toolbox  1.0
A microgeodynamics-based toolkit for rock physics.
 All Classes Files Functions Variables
microgeodynamics.f90
Go to the documentation of this file.
1 
24  USE global
25  TYPE fit
26  REAL(sp)::p1,p2,p3,p4,p5,p6
27  END TYPE fit
28  TYPE melt
29  REAL(sp)::K0,Kp,rho0,dihedral
30  END TYPE melt
32  INTEGER::potential_temperature,basalt_fraction
33  END TYPE composition
34 TYPE mantle
35  TYPE(composition)::comp
36  REAL(sp),DIMENSION(352)::P,Z,rho,vs,vp,K,G,nu,T,X
37  END TYPE mantle
38  TYPE unit_cell
39  REAL(sp)::vs_sol,vp_sol,K_sol,G_sol,nu,Kl,rhol,rho,pressure
40  REAL(sp)::vs_eff,vp_eff,K_eff,G_eff
41  REAL(sp)::vs_obs,vp_obs,K_obs,G_obs
42  REAL(sp)::melt_fraction,contiguity,temperature
43  END TYPE unit_cell
44 
45 CONTAINS
46 
47  FUNCTION clear_unit_cell(cell)
48  IMPLICIT NONE
49  TYPE(unit_cell),INTENT(in)::cell
51 
54  clear_unit_cell%vs_sol=0.0_sp
55  clear_unit_cell%vp_sol=0.0_sp
56  clear_unit_cell%K_sol=0.0_sp
57  clear_unit_cell%G_sol=0.0_sp
58  clear_unit_cell%nu=0.0_sp
59  clear_unit_cell%Kl=0.0_sp
60  clear_unit_cell%rhol=0.0_sp
61  clear_unit_cell%rho=0.0_sp
62  clear_unit_cell%pressure=0.0_sp
63 
64  clear_unit_cell%vs_obs=0.0_sp
65  clear_unit_cell%vp_obs=0.0_sp
66  clear_unit_cell%K_obs=0.0_sp
67  clear_unit_cell%G_obs=0.0_sp
68 
69  clear_unit_cell%vs_eff=0.0_sp
70  clear_unit_cell%vp_eff=0.0_sp
71  clear_unit_cell%K_eff=0.0_sp
72  clear_unit_cell%G_eff=0.0_sp
73 
74  clear_unit_cell%melt_fraction=0.0_sp
75  clear_unit_cell%contiguity=0.0_sp
76  clear_unit_cell%temperature=0.0_sp
77 
78  END FUNCTION clear_unit_cell
79 
80  FUNCTION set_melt(i)
81  IMPLICIT NONE
82  TYPE(melt)::set_melt,m
83  INTEGER,INTENT(in)::i
84 
92 
93  IF(i==2) THEN
94  m%K0=15.5e9_sp
95  m%Kp=7.2_sp
96  m%rho0=2590.0_sp
97  m%dihedral=30.0_sp
98  elseif(i==3) THEN
99  m%K0=24.9e9_sp
100  m%Kp=6.4_sp
101  m%rho0=2610.0_sp
102  m%dihedral=30.0_sp
103  elseif(i==4) THEN
104  m%K0=18.1e9_sp
105  m%Kp=5.5_sp
106  m%rho0=2590.0_sp
107  m%dihedral=30.0_sp
108  ELSE
109  m%K0=16.5e9_sp
110  m%Kp=7.2_sp
111  m%rho0=2610.0_sp
112  m%dihedral=30.0_sp
113  END IF
114  set_melt%K0=m%K0
115  set_melt%kp=m%Kp
116  set_melt%rho0=m%rho0
117  set_melt%dihedral= m%dihedral
118 
119  END FUNCTION set_melt
120 
121  FUNCTION vbw(dihedral,meltfrac)
122  IMPLICIT NONE
123  REAL(sp),INTENT(in)::dihedral,meltfrac
124  REAL(sp)::vbw
125  REAL(sp),DIMENSION(3)::bss,pss,bsl,psl
126  REAL(sp)::b,p,agg,agm
127 
138  bss=(/8.16_sp, -7.71e-2_sp, 1.03e-3_sp/)
139  pss=(/0.424_sp, 9.95e-4_sp, 8.6645e-6_sp/)
140  bsl=(/12.86_sp, -7.85e-2_sp, 1.0043e-3_sp/)
141  psl =(/0.43_sp, 8.63e-5_sp, 2.41e-5_sp/)
142 
143  b=bss(3)*dihedral**2+bss(2)*dihedral+bss(1)
144  p=pss(3)*dihedral**2+pss(2)*dihedral+pss(1)
145 
146  agg=pi-b*meltfrac**p
147 
148  b=bsl(3)*dihedral**2+bsl(2)*dihedral+bsl(1)
149  p=psl(3)*dihedral**2+psl(2)*dihedral+psl(1)
150 
151  agm=b*meltfrac**p
152 
153  vbw=2.0_sp*agg/(2.0_sp*agg+agm)
154  END FUNCTION vbw
155 
156  FUNCTION hmrb(dihedral,meltfrac)
157  IMPLICIT NONE
158  REAL(sp),INTENT(in)::dihedral,meltfrac
159  REAL(sp)::hmrb,temp1,chi1,chi2,agm,agg,temp2
160 
164  !
165  temp2=dihedral*pi/90.0_sp
166  temp1=sqrt(cos(temp2)**2-sin(temp2)*cos(temp2)-0.25_sp*pi+temp2)
167 
168  chi1=abs(0.5_sp*pi -2.0_sp*temp2)/temp1
169 
170  chi2=abs(cos(temp2)-sin(temp2))/temp1
171  agm=chi1*sqrt(meltfrac)
172  agg= 1.0_sp-chi2*sqrt(meltfrac)
173  hmrb=agg/(agg+agm)
174 
175  END FUNCTION hmrb
176 
177  FUNCTION whm(dihedral,meltfrac)
178  IMPLICIT NONE
179  REAL(sp),INTENT(in)::dihedral,meltfrac
180  REAL(sp)::whm
181  REAL(sp),DIMENSION(6)::p
182 
188  p(1)=-8065.0_sp
189  p(2)=6149.0_sp
190  p(3)=-1778.0_sp
191  p(4)=249.0_sp
192  p(5)=-19.77_sp
193  p(6)=1.0_sp
194 
195  whm=p(1)*meltfrac**5+p(2)*meltfrac**4+p(3)*meltfrac**3+p(4)*meltfrac**2+p(5)*meltfrac+p(6)
196  END FUNCTION whm
197 
198  SUBROUTINE vel2mod(vs,vp,rho,G,K,nu)
199  IMPLICIT NONE
200  REAL(sp),INTENT(in)::vs,vp,rho
201  REAL(sp),INTENT(out)::g,k,nu
202 
207 
208  g = rho*vs**2
209  k = rho*(vp**2-1.33_sp*vs**2)
210  nu=(3.0_sp*k-2.0_sp*g)/(6.0_sp*k+2.0_sp*g)
211  END SUBROUTINE vel2mod
212 
213  SUBROUTINE mod2vel(G,K,rho,vs,vp)
214  IMPLICIT NONE
215  REAL(sp),INTENT(in)::g,k,rho
216  REAL(sp),INTENT(out)::vs,vp
217 
219  vs=sqrt(g/rho)
220  vp=sqrt(k/rho+1.33_sp*g/rho)
221 
222  END SUBROUTINE mod2vel
223 
224  FUNCTION elastic_melt(contiguity,meltfrac,K,G,nu,rho,Kl,rhol)
225  REAL(sp),INTENT(in)::contiguity,meltfrac,k,g,nu,kl,rho,rhol
226  REAL(sp),DIMENSION(4)::elastic_melt
227  REAL(sp)::vs1,vp1,k1,g1
228  REAL(sp)::a1,a2,a3,b1,b2,b3
229  REAL(sp)::m,n,h,gg,kboverk,bet,rhobaroverrho,rr0
230 
239 
240  a1=1.8625_sp+0.52594_sp*nu-4.8397_sp*nu**2
241  a2=4.5001_sp-6.1551_sp*nu-4.3634_sp*nu**2
242  a3=-5.6512_sp+6.9159_sp*nu+29.595_sp*nu**2-58.96_sp*nu**3
243  b1=1.6122_sp+0.13527_sp*nu
244  b2=4.5869_sp+3.6086_sp*nu
245  b3=-7.5395_sp-4.8676_sp*nu-4.3182_sp*nu**2
246 
247  m=a1*contiguity+a2*(1.0_sp-contiguity)+a3*(1.0_sp-contiguity)**1.5_sp*contiguity !Equation A5
248  n=b1*contiguity+b2*(1.0_sp-contiguity)+b3*contiguity*(1.0_sp-contiguity)**2 !Equation A6
249  h=1.0_sp-(1.0_sp-contiguity)**m
250  gg=1.0_sp-(1.0_sp-contiguity)**n
251 
252  !normalized bulk modulus of the skeletal framework eq 27, H-m 2008
253  kboverk=(1.0_sp-meltfrac)*h
254  g1=(1.0_sp-meltfrac)*gg
255 
256  bet=k/kl
257  gam=g/k
258  rr0=rhol/rho
259 
260  k1=kboverk+((1.0_sp-kboverk)**2)/(1.0_sp-meltfrac-kboverk+meltfrac*bet)
261  rhobaroverrho=1.0_sp-meltfrac+rr0*meltfrac
262 
263  vs1=sqrt(g1)/sqrt(rhobaroverrho)
264 
265  vp1=sqrt(keffoverk+4.0_sp*gam*g1/3.0_sp)&
266  &/sqrt(1.0_sp+4.0_sp*gam/3.0_sp)/sqrt(rhobaroverrho)
267 
268  elastic_melt(1)=vs1
269  elastic_melt(2)=vp1
270  elastic_melt(3)=k1
271  elastic_melt(4)=g1
272  END FUNCTION elastic_melt
273 
274  SUBROUTINE vinet(K0,rho,Kp,rho0,P,K)
275  IMPLICIT NONE
276  REAL(sp),INTENT(in)::k0,rho,kp,rho0
277  REAL(sp),INTENT(out)::p,k
278  REAL(sp)::z,term1,term2,term3,dpdz
279 
283  z=rho/rho0
284  p=3*k0*z**(2.0_sp/3.0_sp)*(1.0_sp-z**(-1.0_sp/3.0_sp))&
285  &*exp(1.5_sp*(kp-1.0_sp)*(1.0_sp-z**(-1.0_sp/3.0_sp)))
286  term1=3.0_sp*k0*exp(1.5_sp*(kp-1.0_sp)*(1.0_sp-z**(-1.0_sp/3.0_sp)))
287  term2=z**(-1.0_sp/3.0_sp)*(2.0_sp-z**(-1.0_sp/3.0_sp))/3.0_sp
288  term3=0.5_sp*z**(-2.0_sp/3.0_sp)*(1.0_sp-z**(-1.0_sp/3.0_sp))*(kp-1.0_sp)
289  dpdz=term1*(term2+term3)
290  k=z*dpdz
291  END SUBROUTINE vinet
292 
293  SUBROUTINE bm3(K0,rho,Kp,rho0,P,K)
294  IMPLICIT NONE
295  REAL(sp),INTENT(in)::k0,rho,kp,rho0
296  REAL(sp),INTENT(out)::p,k
297  REAL(sp)::z,term1,term2,term3,dpdz
298 
303  z=rho/rho0
304  p=1.5_sp*k0*((rho/rho0)**(7.0_sp/3.0_sp)&
305  &-(rho/rho0)**(5.0_sp/3.0_sp))*&
306  &(1.0_sp+0.75_sp*(kp-4.0_sp)*&
307  &((rho/rho0)**(2.0_sp/3.0_sp)-1.0_sp))
308  term1=0.75_sp*k0*(kp-4.0_sp)*(z**2-z**(4.0_sp/3.0_sp))
309  term2=1.0_sp+0.75_sp*(kp-4.0_sp)*(z**(2.0_sp/3.0_sp)-1.0_sp)
310  term3=0.5_sp*k0*(7.0_sp*z**(4.0_sp/3.0_sp)-5.0_sp*z**(2.0_sp/3.0_sp))
311  dpdz=term1+term2*term3
312  k=z*dpdz
313  END SUBROUTINE bm3
314 
315  FUNCTION load_composition(comp)
316  IMPLICIT NONE
317  TYPE(composition),INTENT(in)::comp
318  TYPE(mantle)::load_composition
319  INTEGER::ii,filesize=352
320  CHARACTER(len=80)::fname
321  REAL(sp)::p,z,t,rho,x,vs,vp,s
322  load_composition%comp=comp
323 !!$ WRITE(fname,'(a15,i2,a1,i4,a4)'),"../input/basalt",comp%basalt_fraction&
324 !!$ &,"_",comp%potential_temperature,".dat"
325  fname='../input/352.dat'
326 
327  OPEN(10,file=fname,form='formatted')
328 
329  DO ii=1,352
330  READ(10,'(8(f10.2,1x))'),p,z,t,rho,x,vs,vp,s
331  load_composition%P(ii)=p*1.0e9_sp! Convert to Pa
332  load_composition%Z(ii)=z*1.0e3_sp! Convert to m
333  load_composition%rho(ii)=rho*1.0e3_sp! Convert to kg/m^3
334  load_composition%vs(ii)=vs*1.0e3_sp! Convert to m/s
335  load_composition%vp(ii)=vp*1.0e3_sp! Convert to m/s
336  load_composition%T(ii)=t !potential temperature
337  load_composition%X(ii)=x !basalt fraction
338  ! Calculate the bulk and shear modulus
339  CALL vel2mod(vs*1.0e3_sp,vp*1.0e3_sp,rho*1.0e3_sp,load_composition%G(ii),&
340  &load_composition%K(ii),load_composition%nu(ii))
341  END DO
342  close(10)
343  END FUNCTION load_composition
344 
345  FUNCTION bisection(x1,x2,eps,cell,magma)
346  IMPLICIT NONE
347  REAL(sp),INTENT(in)::x1,x2,eps
348  TYPE(unit_cell),INTENT(in)::cell
349  TYPE(melt),INTENT(in)::magma
350  REAL(sp)::bisection,dx,f,fmid,xmid
351  REAL(sp),DIMENSION(2)::temp
352  INTEGER::ii,maxit,jj
353  maxit=200
354  temp=funcd(x2,cell,magma)
355  fmid=temp(1)
356  temp=funcd(x1,cell,magma)
357  f=temp(1)
358  IF(f*fmid>0) THEN
359  print*,'root not bracketted in bisection method'
360  bisection=x1
361  return
362  END IF
363  IF (f<0.0_sp) THEN
364  bisection=x1
365  dx=x2-x1
366  ELSE
367  bisection=x2
368  dx=x1-x2
369  END IF
370  DO jj=1,maxit
371  dx=0.5_sp*dx
372  xmid=bisection+dx
373  temp=funcd(xmid,cell,magma)
374  fmid=temp(1)
375  IF (fmid<0.0_sp) THEN
376  bisection=xmid
377  END IF
378 
379  IF(abs(dx)<eps.OR.fmid==0.0_sp) return
380 
381  END DO
382 
383  END FUNCTION bisection
384 
385  FUNCTION nonlinear_solve(x1,x2,eps,cell,magma)
386  IMPLICIT NONE
387  REAL(sp),INTENT(in)::x1,x2,eps
388  TYPE(unit_cell),INTENT(in)::cell
389  TYPE(melt),INTENT(in)::magma
390  REAL(sp)::sol,nonlinear_solve,xlo,xhi,dx,dxold,temp
391  REAL(sp),DIMENSION(2)::flo,fhi,f
392  INTEGER::ii,maxit,jj
393 
400  flo=funcd(x1,cell,magma)
401  fhi=funcd(x2,cell,magma)
402 !!$
403 
404  !! Does the function cross zero within the provided bounds?
405  IF((flo(1)<0.0_sp.AND.fhi(1)<0.0_sp).OR.(flo(1)>0.0_sp.AND.fhi(1)>0.0_sp)) THEN
406  print*,'nonlinear solve failed, function not bracketted within bounds'
407  return
408  elseif(flo(1)==0.0_sp) THEN
409  nonlinear_solve=x1 ! The lower limit is the solution
410  elseif(fhi(1)==0.0_sp) THEN
411  nonlinear_solve=x2 ! The upper limit is the solution
412  elseif(flo(1)<0.0_sp) THEN
413  xlo=x1 ! Set it such that f(xlo)<0
414  xhi=x2
415  ELSE
416  xhi=x1
417  xlo=x2
418  END IF
419  nonlinear_solve=0.5_sp*(x1+x2) ! Start with bisection
420  dxold=abs(x2-x1) !stepsize before last stepsize
421  dx=dxold
422  f=funcd(nonlinear_solve,cell,magma)! Evaluate function at the current location
423 
424  maxit=100
425  DO jj=1,maxit
426  !! Bisect if the Newton iteration is
427  !! out of range or not decreasing fast enough
428  IF(((nonlinear_solve-xhi)*f(2)-f(1))*((nonlinear_solve-xlo)*f(2)-f(1))>0.0_sp.OR.&
429  &abs(2.0*f(1))>abs(dxold*f(2))) THEN
430  dxold=dx
431  dx=0.5_sp*(xhi-xlo)
432  nonlinear_solve=xlo+dx
433 
434  IF(xlo==nonlinear_solve) return
435  ELSE
436  dxold=dx
437  dx=f(1)/f(2)
438  temp=nonlinear_solve
440 
441  IF(temp==nonlinear_solve) return
442  END IF
443  IF(abs(dx)<eps) return
444 
445  f=funcd(nonlinear_solve,cell,magma)
446  IF(f(1)<0.0_sp) THEN
447  xlo=nonlinear_solve
448  ELSE
449  xhi=nonlinear_solve
450  END IF
451  END DO
452  END FUNCTION nonlinear_solve
453 
454 
455  FUNCTION funcd(x,cell,magma)
456  IMPLICIT NONE
457  REAL(sp),INTENT(in)::x
458  TYPE(unit_cell),INTENT(in)::cell
459  TYPE(melt),INTENT(in)::magma
460  REAL(sp),DIMENSION(2)::funcd
461  REAL(sp)::dx,cont
462  REAL(sp),DIMENSION(4)::effective
463 
468 
469  cont=whm(magma%dihedral,x)
470  effective= elastic_melt(cont,x,cell%K_sol,cell%G_sol&
471  &,cell%nu,cell%rho,cell%Kl,cell%rhol)
472 
473  funcd(1)=effective(1)*cell%vs_sol-cell%vs_obs
474  dx=5.0e-3_sp
475  cont=whm(magma%dihedral,x+dx)
476  effective= elastic_melt(cont,x+dx,cell%K_sol,cell%G_sol&
477  &,cell%nu,cell%rho,cell%Kl,cell%rhol)
478  funcd(2)=(effective(1)*cell%vs_sol-cell%vs_obs-funcd(1))/dx
479 
480  END FUNCTION funcd
481 
482  FUNCTION temp2vs(T,basalt_fraction)!<model)
483  REAL(sp),INTENT(in)::t
484  INTEGER,INTENT(in)::basalt_fraction
485  !type(mantle),intent(in)::model
486  TYPE(fit),DIMENSION(3)::f
487  REAL(sp)::temp2vs,temptol
488  integer::ii,n
489 
491  f=temperaturefit(basalt_fraction)
492  temp2vs=f(2)%p1*t**2+f(2)%p2*t+f(2)%p3
493  temp2vs=temp2vs*1.0e3_sp
494  END FUNCTION temp2vs
495 
496  FUNCTION temp2vp(T,basalt_fraction)
497 
499  REAL(sp),INTENT(in)::t
500  INTEGER,INTENT(in)::basalt_fraction
501  TYPE(fit),DIMENSION(3)::f
502  REAL(sp)::temp2vp
503  f=temperaturefit(basalt_fraction)
504  temp2vp=f(3)%p1*t**2+f(3)%p2*t+f(3)%p3
505  temp2vp=temp2vp*1.0e3_sp
506  END FUNCTION temp2vp
507 
508  FUNCTION temp2rho(T,basalt_fraction)
509 
511  REAL(sp),INTENT(in)::t
512  INTEGER,INTENT(in)::basalt_fraction
513  TYPE(fit),DIMENSION(3)::f
514  REAL(sp)::temp2rho
515  f=temperaturefit(basalt_fraction)
516  temp2rho=f(1)%p1*t**2+f(1)%p2*t+f(1)%p3
517  temp2rho=temp2rho*1.0e3_sp
518  END FUNCTION temp2rho
519 
520  FUNCTION temperaturefit(basalt_fraction)
521  INTEGER,INTENT(in)::basalt_fraction
522  TYPE(fit),DIMENSION(3)::f,temperaturefit
523  integer::ii,n
524  real(sp),dimension(22,10)::par
525  real(sp)::x
526 
532 
533 
534  ! Enter coefficients for polynomial fit for vs
535  open(10,file='../input/xu_fit_param.csv',form='formatted')
536  x=dble(basalt_fraction)/100.0_sp
537  n=22
538  do ii=1,22
539  read(10,'(10(1x,es10.3))'),par(ii,:)
540  end do
541  close(10)
542  do ii=1,22
543  if(x==par(ii,1)) then
544  !Density
545  f(1)%p1=par(ii,2)
546  f(1)%p2=par(ii,3)
547  f(1)%p3=par(ii,4)
548  !Vs
549  f(2)%p1=par(ii,5)
550  f(2)%p2=par(ii,6)
551  f(2)%p3=par(ii,7)
552  !Vp
553  f(3)%p1=par(ii,8)
554  f(3)%p2=par(ii,9)
555  f(3)%p3=par(ii,10)
556  end if
557  end do
559  END FUNCTION temperaturefit
560 
561  FUNCTION depth2dt(depth)
562  IMPLICIT NONE
563  REAL(sp),INTENT(in)::depth
564  REAL(sp)::depth2dt,dzdt
565 
568  dzdt=0.1_sp ! km/K
569  depth2dt=(depth-410.0_sp)/dzdt
570  END FUNCTION depth2dt
571 
572 END MODULE microgeodynamics