endgame_Fortran/functions.F90 at master · arnaudvericel/endgame_Fortran · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
!- Mathematical functions used by endgame
module functions
 use config,    only:r0,T0,cs0,sigma_g0,p,q,mstar,rho_g0,eta0,rbump,epsi0,ibump,phi,w,epsimax,alpha
 use config,     only:gg,pi,au,earthr,ibr

 implicit none

 public         :: sigma_g,Temp,omega_k,cs,vdrift,vvisc,rho_g,vk,hoverr,h,epsi

 contains

real function sigma_g(r)
 real, intent(in) :: r

 sigma_g = sigma_g0*(r/r0)**(-p) + ibump*phi*sigma_g0*exp(-(r-rbump)**2/(2*w**2))
 return
end function sigma_g

real function rho_g(r)
 real, intent(in) :: r

 rho_g = sigma_g(r)/(sqrt(2*pi)*h(r))
 return
end function rho_g

real function Temp(r)
 real, intent(in) :: r

 Temp = T0*(r/r0)**(-q)
 return
end function Temp

real function omega_k(r)
 real, intent(in) :: r

 omega_k = sqrt(gg*mstar/r**3)
 return
end function omega_k

real function cs(r)
 real, intent(in) :: r

 cs = cs0*(r/r0)**(-q/2.)
 return
end function cs

real function vk(r)
 real, intent(in) :: r

 vk = omega_k(r)*r
 return
end function vk

real function vdrift(St,r)
 real, intent(in) :: St,r
 real             :: deriv,eps

 deriv = hoverr(r)**2 * vk(r) * ((press(r+earthr) - press(r-earthr))/(2*earthr)) * r/press(r)

 if (ibr==1) then
    eps=epsi(r)
 else
    eps = 0.
 endif

 vdrift = St / ((1+eps)**2 + St**2) * deriv
 return
end function vdrift

real function vvisc(St,r)
 real, intent(in) :: St,r
 real             :: K, deriv,eps

 K = 3/(r*rho_g(r)*vk(r))
 deriv = (rho_g(r+earthr)*nu(r+earthr)*(r+earthr)*vk(r+earthr) - rho_g(r-earthr)*nu(r-earthr)*(r-earthr)*vk(r-earthr)) / (2*earthr)

 if (ibr==1) then
    eps=epsi(r)
 else
    eps=0.
 endif

 vvisc = (1+eps)/((1+eps)**2 + St**2) * K * deriv
 return
end function vvisc

real function nu(r)
 real, intent(in) :: r
 nu = alpha*cs(r)*h(r)
 return
end function nu

real function hoverr(r)
 real, intent(in) :: r

 hoverr = h(r)/r
 return
end function hoverr

real function h(r)
 real, intent(in) :: r

 h = cs(r)/omega_k(r)
 return
end function h

real function press(r)
 real, intent(in) :: r

 press = cs(r)**2*rho_g(r)
 return
end function press

real function epsi(r)
 real, intent(in) :: r

 epsi = epsi0 + ibump*(epsimax-epsi0)*exp(-(r-rbump)**2/(2*w**2))
 return
end function epsi

end module functions