-
Notifications
You must be signed in to change notification settings - Fork 14
/
distn.f
95 lines (90 loc) · 3.1 KB
/
distn.f
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
c
c
c
subroutine distn(nv,nt,v0,thet0,ff,dfdx,dfdy,denn,den,ird,cdvt2)
implicit double precision (a-h,o-z)
c include 'clich1.h'
dimension den(ird)
c Could speed up by eliminating repeated searches of grid
c in terp2, for x and y.
c
c interpolate (using density) to get distribution functions
c Location alogrithm for interpolation points assumes that density
c is a decreasing function of radius.
c
do 1 i=1,ird
ir=i
if(denn.gt.den(ir)) goto 2
1 continue
2 idm=nva
c
c
ff=(denn/den(ir))*terp2p(v0,thet0,nv,vs,nt,thets,
1 f(1,1,ir),fxx(1,1,ir),fyy(1,1,ir),fxxyy(1,1,ir),idm,0,0,1)
dfdx=(denn/den(ir))*terp2p(v0,thet0,nv,vs,nt,thets,
1 f(1,1,ir),fxx(1,1,ir),fyy(1,1,ir),fxxyy(1,1,ir),idm,1,0,0)
dfdy=(denn/den(ir))*terp2p(v0,thet0,nv,vs,nt,thets,
1 f(1,1,ir),fxx(1,1,ir),fyy(1,1,ir),fxxyy(1,1,ir),idm,0,1,0)
c
if(ir.eq.1.or.(ir.eq.ird.and.denn.le.den(ird))) return
im=ir-1
ffm=(denn/den(im))*terp2p(v0,thet0,nv,vs,nt,thets,
1 f(1,1,im),fxx(1,1,im),fyy(1,1,im),fxxyy(1,1,im),idm,0,0,0)
dfdxm=(denn/den(im))*terp2p(v0,thet0,nv,vs,nt,thets,
1 f(1,1,im),fxx(1,1,im),fyy(1,1,im),fxxyy(1,1,im),idm,1,0,0)
dfdym=(denn/den(im))*terp2p(v0,thet0,nv,vs,nt,thets,
1 f(1,1,im),fxx(1,1,im),fyy(1,1,im),fxxyy(1,1,im),idm,0,1,0)
c
ddm=denn-den(im)
ddp=den(ir)-denn
dd=ddm+ddp
ff=( ffm*ddp + ff*ddm ) / dd
dfdx=( dfdxm*ddp + dfdx*ddm ) / dd
dfdy=( dfdym*ddp + dfdy*ddm ) / dd
c
return
end
c!
subroutine distr(nv,nt,v0,thet0,ff,dfdx,dfdy,rho,rovera,ird,cdvt2)
implicit double precision (a-h,o-z)
include 'clich1.h'
dimension den(ird)
c Could speed up by eliminating repeated searches of grid
c in terp2, for x and y.
c
c interpolate (using density) to get distribution functions
c Location alogrithm for interpolation points assumes that density
c is a decreasing function of radius.
c
do 1 i=1,ird
ir=i
if(denn.gt.den(ir)) goto 2
1 continue
2 idm=nva
c
c
ff=(denn/den(ir))*terp2p(v0,thet0,nv,vs,nt,thets,
1 f(1,1,ir),fxx(1,1,ir),fyy(1,1,ir),fxxyy(1,1,ir),idm,0,0,1)
dfdx=(denn/den(ir))*terp2p(v0,thet0,nv,vs,nt,thets,
1 f(1,1,ir),fxx(1,1,ir),fyy(1,1,ir),fxxyy(1,1,ir),idm,1,0,0)
dfdy=(denn/den(ir))*terp2p(v0,thet0,nv,vs,nt,thets,
1 f(1,1,ir),fxx(1,1,ir),fyy(1,1,ir),fxxyy(1,1,ir),idm,0,1,0)
c
if(ir.eq.1.or.(ir.eq.ird.and.denn.le.den(ird))) return
im=ir-1
ffm=(denn/den(im))*terp2p(v0,thet0,nv,vs,nt,thets,
1 f(1,1,im),fxx(1,1,im),fyy(1,1,im),fxxyy(1,1,im),idm,0,0,0)
dfdxm=(denn/den(im))*terp2p(v0,thet0,nv,vs,nt,thets,
1 f(1,1,im),fxx(1,1,im),fyy(1,1,im),fxxyy(1,1,im),idm,1,0,0)
dfdym=(denn/den(im))*terp2p(v0,thet0,nv,vs,nt,thets,
1 f(1,1,im),fxx(1,1,im),fyy(1,1,im),fxxyy(1,1,im),idm,0,1,0)
c
ddm=denn-den(im)
ddp=den(ir)-denn
dd=ddm+ddp
ff=( ffm*ddp + ff*ddm ) / dd
dfdx=( dfdxm*ddp + dfdx*ddm ) / dd
dfdy=( dfdym*ddp + dfdy*ddm ) / dd
c
return
end