-
Notifications
You must be signed in to change notification settings - Fork 2
/
adams.f90
86 lines (68 loc) · 1.46 KB
/
adams.f90
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
subroutine adams
! Adams-Bashforth scheme
use vars
use params, only: dowallx, dowally
implicit none
real dtdx, dtdy, dtdz, rhox, rhoy, rhoz , a1, a2
integer i,j,k
call t_startf('adams')
u1(:,:,:) = u(:,:,:)
v1(:,:,:) = v(:,:,:)
w1(:,:,:) = w(:,:,:)
do k=1,nzm
do j=1,ny
do i=1,nx
u(i,j,k) = u1(i,j,k)+dt3(na) &
*(at*dudt(i,j,k,na)+bt*dudt(i,j,k,nb)+ct*dudt(i,j,k,nc))
v(i,j,k) = v1(i,j,k)+dt3(na) &
*(at*dvdt(i,j,k,na)+bt*dvdt(i,j,k,nb)+ct*dvdt(i,j,k,nc))
w(i,j,k) = w1(i,j,k)+dt3(na) &
*(at*dwdt(i,j,k,na)+bt*dwdt(i,j,k,nb)+ct*dwdt(i,j,k,nc))
end do
end do
end do
! Exchange boundaries:
call boundaries(1)
! compute time averaged velocties for second-order advection of scalars:
dtdx = dtn/dx
dtdy = dtn/dy
dtdz = dtn/dz
a1 = 0.5
a2 = 0.5
if(nstep.eq.1) then
a1 = 1.
a2 = 0.
end if
do k=1,nzm
rhox = rho(k)*dtdx
do j=dimy1_u,dimy2_u
do i=dimx1_u,dimx2_u
u1(i,j,k) = (a1*u(i,j,k)+a2*u1(i,j,k))*rhox
end do
end do
end do
do k=1,nzm
rhoy = rho(k)*dtdy
do j=dimy1_v,dimy2_v
do i=dimx1_v,dimx2_v
v1(i,j,k) = (a1*v(i,j,k)+a2*v1(i,j,k))*rhoy
end do
end do
end do
do k=1,nzm
do j=1,ny
do i=1,nx
misc(i,j,k) = (a1*w(i,j,k)+a2*w1(i,j,k))
end do
end do
end do
do k=1,nzm
rhoz = rhow(k)*dtdz
do j=dimy1_w,dimy2_w
do i=dimx1_w,dimx2_w
w1(i,j,k) = (a1*w(i,j,k)+a2*w1(i,j,k))*rhoz
end do
end do
end do
call t_stopf ('adams')
end subroutine adams