source: palm/trunk/SOURCE/init_rankine.f90 @ 107

Last change on this file since 107 was 107, checked in by raasch, 14 years ago

further bugfix for non-cyclic BCs

  • Property svn:keywords set to Id
File size: 4.0 KB
Line 
1 SUBROUTINE init_rankine
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6! Initial profiles are reset to constant profiles
7!
8! Former revisions:
9! -----------------
10! $Id: init_rankine.f90 107 2007-08-17 13:54:45Z raasch $
11!
12! 75 2007-03-22 09:54:05Z raasch
13! uxrp, vynp eliminated, 2nd+3rd argument removed from exchange horiz
14!
15! RCS Log replace by Id keyword, revision history cleaned up
16!
17! Revision 1.11  2005/03/26 20:38:49  raasch
18! Arguments for non-cyclic boundary conditions added to argument list of
19! routine exchange_horiz
20!
21! Revision 1.1  1997/08/11 06:18:43  raasch
22! Initial revision
23!
24!
25! Description:
26! ------------
27! Initialize a (nondivergent) Rankine eddy with a vertical axis in order to test
28! the advection terms and the pressure solver.
29!------------------------------------------------------------------------------!
30
31    USE arrays_3d
32    USE constants
33    USE grid_variables
34    USE indices
35    USE control_parameters
36
37    IMPLICIT NONE
38
39    INTEGER ::  i, ic, j, jc, k, kc1, kc2
40    REAL    ::  alpha, betrag, radius, rc, uw, vw, x, y
41
42!
43!-- Default: eddy radius rc, eddy strength z,
44!--          position of eddy centre: ic, jc, kc1, kc2
45    rc  =  4.0 * dx
46    ic  =  ( nx+1 ) / 2
47    jc  =  ic
48    kc1 = nzb
49    kc2 = nzt+1
50
51!
52!-- Reset initial profiles to constant profiles
53    IF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 )  THEN
54       DO  i = nxl-1, nxr+1
55          DO  j = nys-1, nyn+1
56             pt(:,j,i) = pt_init
57             u(:,j,i)  = u_init
58             v(:,j,i)  = v_init
59          ENDDO
60       ENDDO
61    ENDIF
62
63!
64!-- Compute the u-component.
65    DO  i = nxl, nxr
66       DO  j = nys, nyn
67          x = ( i - ic - 0.5 ) * dx
68          y = ( j - jc ) * dy
69          radius = SQRT( x**2 + y**2 )
70          IF ( radius <= 2.0 * rc )  THEN
71             betrag = radius / ( 2.0 * rc ) * 0.08
72          ELSEIF ( radius > 2.0 * rc  .AND.  radius < 8.0 * rc )  THEN
73             betrag = 0.08 * EXP( -( radius - 2.0 * rc ) / 2.0 )
74          ELSE
75             betrag = 0.0
76          ENDIF
77
78          IF ( x == 0.0 )  THEN
79             IF ( y > 0.0 )  THEN
80                alpha = pi / 2.0
81             ELSEIF ( y < 0.0 )  THEN
82                alpha = 3.0 * pi / 2.0
83             ENDIF
84          ELSE
85             IF ( x < 0.0 )  THEN
86                alpha = ATAN( y / x ) + pi
87             ELSE
88                IF ( y < 0.0 )  THEN
89                   alpha = ATAN( y / x ) + 2.0 * pi
90                ELSE
91                   alpha = ATAN( y / x )
92                ENDIF
93             ENDIF
94          ENDIF
95
96          uw = -SIN( alpha ) * betrag
97
98          DO  k = kc1, kc2
99             u(k,j,i) = u(k,j,i) + uw
100          ENDDO
101       ENDDO
102    ENDDO
103
104!
105!-- Compute the v-component.
106    DO  i = nxl, nxr
107       DO  j = nys, nyn
108          x = ( i - ic ) * dx
109          y = ( j - jc - 0.5) * dy
110          radius = SQRT( x**2 + y**2 )
111          IF ( radius <= 2.0 * rc )  THEN
112             betrag = radius / ( 2.0 * rc ) * 0.08
113          ELSEIF ( radius > 2.0 * rc  .AND.  radius < 8.0 * rc )  THEN
114             betrag = 0.08 * EXP( -( radius - 2.0 * rc ) / 2.0 )
115          ELSE
116             betrag = 0.0
117          ENDIF
118
119          IF ( x == 0.0 )  THEN
120             IF ( y > 0.0 )  THEN
121                alpha = pi / 2.0
122             ELSEIF ( y < 0.0 )  THEN
123                alpha = 3.0 * pi / 2.0
124             ENDIF
125          ELSE
126             IF ( x < 0.0 )  THEN
127                alpha = ATAN( y / x ) + pi
128             ELSE
129                IF ( y < 0.0 )  THEN
130                   alpha = ATAN( y / x ) + 2.0 * pi
131                ELSE
132                   alpha = ATAN( y / x )
133                ENDIF
134             ENDIF
135          ENDIF
136
137          vw = COS( alpha ) * betrag
138
139          DO  k = kc1, kc2
140             v(k,j,i) = v(k,j,i) + vw
141          ENDDO
142       ENDDO
143    ENDDO
144
145!
146!-- Exchange of boundary values for the velocities.
147    CALL exchange_horiz( u )
148    CALL exchange_horiz( v )
149!
150!-- Make velocity field nondivergent.
151    n_sor = nsor_ini
152    CALL pres
153    n_sor = nsor
154
155 END SUBROUTINE init_rankine
Note: See TracBrowser for help on using the repository browser.