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

Last change on this file since 1 was 1, checked in by raasch, 15 years ago

Initial repository layout and content

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