source: palm/trunk/SOURCE/init_advec.f90 @ 347

Last change on this file since 347 was 4, checked in by raasch, 18 years ago

Id keyword set as property for all *.f90 files

  • Property svn:keywords set to Id
File size: 6.8 KB
Line 
1 SUBROUTINE init_advec
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: init_advec.f90 4 2007-02-13 11:33:16Z heinze $
11! RCS Log replace by Id keyword, revision history cleaned up
12!
13! Revision 1.6  2004/04/30 11:59:31  raasch
14! impulse_advec renamed momentum_advec
15!
16! Revision 1.1  1999/02/05 09:07:38  raasch
17! Initial revision
18!
19!
20! Description:
21! ------------
22! Initialize constant coefficients and parameters for certain advection schemes.
23!------------------------------------------------------------------------------!
24
25    USE advection
26    USE arrays_3d
27    USE indices
28    USE control_parameters
29
30    IMPLICIT NONE
31
32    INTEGER :: i, intervals, j, k
33    REAL    :: delt, dn, dnneu, ex1, ex2, ex3, ex4, ex5, ex6, spl_alpha, &
34               spl_beta, sterm
35    REAL, DIMENSION(:), ALLOCATABLE ::  spl_u, temp
36
37
38    IF ( scalar_advec == 'bc-scheme' )  THEN
39
40!
41!--    Compute exponential coefficients for the Bott-Chlond scheme
42       intervals = 1000
43       ALLOCATE( aex(intervals), bex(intervals), dex(intervals), eex(intervals) )
44
45       delt  = 1.0 / REAL( intervals )
46       sterm = delt * 0.5
47
48       DO  i = 1, intervals
49
50          IF ( sterm > 0.5 )  THEN
51             dn = -5.0
52          ELSE
53             dn = 5.0
54          ENDIF
55
56          DO  j = 1, 15
57             ex1 = dn * EXP( -dn ) - EXP( 0.5 * dn ) + EXP( -0.5 * dn )
58             ex2 = EXP( dn ) - EXP( -dn )
59             ex3 = EXP( -dn ) * ( 1.0 - dn ) - 0.5 * EXP(  0.5 * dn ) &
60                                             - 0.5 * EXP( -0.5 * dn )
61             ex4 = EXP( dn ) + EXP( -dn )
62             ex5 = dn * sterm + ex1 / ex2
63             ex6 = sterm + ( ex3 * ex2 - ex4 * ex1 ) / ( ex2 * ex2 )
64             dnneu = dn - ex5 / ex6
65             dn  = dnneu
66          ENDDO
67
68          IF ( sterm < 0.5 )  dn = MAX(  2.95E-2, dn )
69          IF ( sterm > 0.5 )  dn = MIN( -2.95E-2, dn )
70          ex1 = EXP( -dn )
71          ex2 = EXP( dn ) - ex1
72          aex(i) = -ex1 / ex2
73          bex(i) = 1.0 / ex2
74          dex(i) = dn
75          eex(i) = EXP( dex(i) * 0.5 )
76          sterm = sterm + delt
77
78       ENDDO
79
80    ENDIF
81
82    IF ( momentum_advec == 'ups-scheme' .OR. scalar_advec  == 'ups-scheme' )  &
83    THEN
84
85!
86!--    Provide the constant parameters for the Upstream-Spline advection scheme.
87!--    In x- und y-direction the Sherman-Morrison formula is applied
88!--    (cf. Press et al, 1986 (Numerical Recipes)).
89!
90!--    Allocate nonlocal arrays
91       ALLOCATE( spl_z_x(0:nx), spl_z_y(0:ny), spl_tri_x(5,0:nx), &
92                 spl_tri_y(5,0:ny), spl_tri_zu(5,nzb:nzt+1),      &
93                 spl_tri_zw(5,nzb:nzt+1) )
94
95!
96!--    Provide diagonal elements of the tridiagonal matrices for all
97!--    directions
98       spl_tri_x(1,:)  = 2.0
99       spl_tri_y(1,:)  = 2.0
100       spl_tri_zu(1,:) = 2.0
101       spl_tri_zw(1,:) = 2.0
102
103!
104!--    Elements of the cyclic tridiagonal matrix
105!--    (same for all horizontal directions)
106       spl_alpha = 0.5
107       spl_beta  = 0.5
108
109!
110!--    Sub- and superdiagonal elements, x-direction
111       spl_tri_x(2,0:nx) = 0.5
112       spl_tri_x(3,0:nx) = 0.5
113
114!
115!--    mMdify the diagonal elements (Sherman-Morrison)
116       spl_gamma_x    = -spl_tri_x(1,0)
117       spl_tri_x(1,0)  = spl_tri_x(1,0) - spl_gamma_x
118       spl_tri_x(1,nx) = spl_tri_x(1,nx) - spl_alpha * spl_beta / spl_gamma_x
119
120!
121!--    Split the tridiagonal matrix for Thomas algorithm
122       spl_tri_x(4,0) = spl_tri_x(1,0)
123       DO  i = 1, nx
124          spl_tri_x(5,i) = spl_tri_x(2,i) / spl_tri_x(4,i-1)
125          spl_tri_x(4,i) = spl_tri_x(1,i) - spl_tri_x(5,i) * spl_tri_x(3,i-1)
126       ENDDO
127
128!
129!--    Allocate arrays required locally
130       ALLOCATE( temp(0:nx), spl_u(0:nx) )
131
132!
133!--    Provide "corrective vector", x-direction
134       spl_u(0)      = spl_gamma_x
135       spl_u(1:nx-1) = 0.0
136       spl_u(nx)     = spl_alpha
137
138!
139!--    Solve the system of equations for the corrective vector
140!--    (Sherman-Morrison)
141       temp(0) = spl_u(0)
142       DO  i = 1, nx
143          temp(i) = spl_u(i) - spl_tri_x(5,i) * temp(i-1)
144       ENDDO
145       spl_z_x(nx) = temp(nx) / spl_tri_x(4,nx)
146       DO  i = nx-1, 0, -1
147          spl_z_x(i) = ( temp(i) - spl_tri_x(3,i) * spl_z_x(i+1) ) / &
148                       spl_tri_x(4,i)
149       ENDDO
150
151!
152!--    Deallocate local arrays, for they are allocated in a different way for
153!--    operations in y-direction
154       DEALLOCATE( temp, spl_u )   
155   
156!
157!--    Provide sub- and superdiagonal elements, y-direction
158       spl_tri_y(2,0:ny) = 0.5
159       spl_tri_y(3,0:ny) = 0.5
160
161!
162!--    Modify the diagonal elements (Sherman-Morrison)
163       spl_gamma_y    = -spl_tri_y(1,0) 
164       spl_tri_y(1,0)  = spl_tri_y(1,0) - spl_gamma_y
165       spl_tri_y(1,ny) = spl_tri_y(1,ny) - spl_alpha * spl_beta / spl_gamma_y
166
167!
168!--    Split the tridiagonal matrix for Thomas algorithm
169       spl_tri_y(4,0) = spl_tri_y(1,0)
170       DO  j = 1, ny
171          spl_tri_y(5,j) = spl_tri_y(2,j) / spl_tri_y(4,j-1)
172          spl_tri_y(4,j) = spl_tri_y(1,j) - spl_tri_y(5,j) * spl_tri_y(3,j-1)
173       ENDDO
174
175!
176!--    Allocate arrays required locally
177       ALLOCATE( temp(0:ny), spl_u(0:ny) )
178
179!
180!--    Provide "corrective vector", y-direction
181       spl_u(0)      = spl_gamma_y
182       spl_u(1:ny-1) = 0.0
183       spl_u(ny)     = spl_alpha
184   
185!
186!--    Solve the system of equations for the corrective vector
187!--    (Sherman-Morrison)
188       temp = 0.0
189       spl_z_y = 0.0
190       temp(0) = spl_u(0)
191       DO  j = 1, ny
192          temp(j) = spl_u(j) - spl_tri_y(5,j) * temp(j-1)
193       ENDDO
194       spl_z_y(ny) = temp(ny) / spl_tri_y(4,ny)
195       DO  j = ny-1, 0, -1
196          spl_z_y(j) = ( temp(j) - spl_tri_y(3,j) * spl_z_y(j+1) ) / &
197                       spl_tri_y(4,j)
198       ENDDO
199
200!
201!--    deallocate local arrays, for they are no longer required
202       DEALLOCATE( temp, spl_u )
203
204!
205!--    provide sub- and superdiagonal elements, z-direction
206       spl_tri_zu(2,nzb)   = 0.0
207       spl_tri_zu(2,nzt+1) = 1.0
208       spl_tri_zw(2,nzb)   = 0.0
209       spl_tri_zw(2,nzt+1) = 1.0
210
211       spl_tri_zu(3,nzb)   = 1.0
212       spl_tri_zu(3,nzt+1) = 0.0
213       spl_tri_zw(3,nzb)   = 1.0
214       spl_tri_zw(3,nzt+1) = 0.0
215
216       DO  k = nzb+1, nzt
217          spl_tri_zu(2,k) = dzu(k) / ( dzu(k) + dzu(k+1) )
218          spl_tri_zw(2,k) = dzw(k) / ( dzw(k) + dzw(k+1) )
219          spl_tri_zu(3,k) = 1.0 - spl_tri_zu(2,k)
220          spl_tri_zw(3,k) = 1.0 - spl_tri_zw(2,k)
221       ENDDO
222   
223       spl_tri_zu(4,nzb) = spl_tri_zu(1,nzb)
224       spl_tri_zw(4,nzb) = spl_tri_zw(1,nzb)
225       DO  k = nzb+1, nzt+1
226          spl_tri_zu(5,k) = spl_tri_zu(2,k) / spl_tri_zu(4,k-1)
227          spl_tri_zw(5,k) = spl_tri_zw(2,k) / spl_tri_zw(4,k-1)
228          spl_tri_zu(4,k) = spl_tri_zu(1,k) - spl_tri_zu(5,k) * &
229                            spl_tri_zu(3,k-1)
230          spl_tri_zw(4,k) = spl_tri_zw(1,k) - spl_tri_zw(5,k) * &
231                            spl_tri_zw(3,k-1)
232       ENDDO
233
234    ENDIF
235
236 END SUBROUTINE init_advec
Note: See TracBrowser for help on using the repository browser.