source: palm/trunk/SOURCE/init_rankine.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: 4.3 KB
Line 
1 SUBROUTINE init_rankine
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Log: init_rankine.f90,v $
11! Revision 1.11  2005/03/26 20:38:49  raasch
12! Arguments for non-cyclic boundary conditions added to argument list of
13! routine exchange_horiz
14!
15! Revision 1.10  2004/01/28 15:11:27  raasch
16! Setting of bt before calling pres removed
17!
18! Revision 1.9  2001/03/30 07:32:49  raasch
19! Arguments removed from pres,
20! Translation of remaining German identifiers (variables, subroutines, etc.)
21!
22! Revision 1.8  2001/01/22 19:44:52  letzel
23! All comments translated into English.
24!
25! Revision 1.7  2001/01/22 07:20:30  raasch
26! Module test_variables removed
27!
28! Revision 1.6  1998/07/06 12:17:44  raasch
29! + USE test_variables
30!
31! Revision 1.5  1998/04/15 11:20:38  raasch
32! pt, usws und vsws werden an pres uebergeben
33!
34! Revision 1.4  1997/09/12 06:26:30  raasch
35! Wirbelradius jetzt proportional zur Gitterweite
36!
37! Revision 1.3  1997/09/09 08:29:24  raasch
38! Kehrwerte der Gitterweiten implementiert
39!
40! Revision 1.2  1997/08/26 06:31:33  raasch
41! +tend in pres-Aufruf
42!
43! Revision 1.1  1997/08/11 06:18:43  raasch
44! Initial revision
45!
46!
47! Description:
48! ------------
49! Initialize a (nondivergent) Rankine eddy with a vertical axis in order to test
50! the advection terms and the pressure solver.
51!------------------------------------------------------------------------------!
52
53    USE arrays_3d
54    USE constants
55    USE grid_variables
56    USE indices
57    USE control_parameters
58
59    IMPLICIT NONE
60
61    INTEGER ::  i, ic, j, jc, k, kc1, kc2
62    REAL    ::  alpha, betrag, radius, rc, uw, vw, x, y
63
64!
65!-- Default: eddy radius rc, eddy strength z,
66!--          position of eddy centre: ic, jc, kc1, kc2
67    rc  =  4.0 * dx
68    ic  =  ( nx+1 ) / 2
69    jc  =  ic
70    kc1 = nzb
71    kc2 = nzt+1
72
73!
74!-- Compute the u-component.
75    DO  i = nxl, nxr
76       DO  j = nys, nyn
77          x = ( i - ic - 0.5 ) * dx
78          y = ( j - jc ) * dy
79          radius = SQRT( x**2 + y**2 )
80          IF ( radius <= 2.0 * rc )  THEN
81             betrag = radius / ( 2.0 * rc ) * 0.08
82          ELSEIF ( radius > 2.0 * rc  .AND.  radius < 8.0 * rc )  THEN
83             betrag = 0.08 * EXP( -( radius - 2.0 * rc ) / 2.0 )
84          ELSE
85             betrag = 0.0
86          ENDIF
87
88          IF ( x == 0.0 )  THEN
89             IF ( y > 0.0 )  THEN
90                alpha = pi / 2.0
91             ELSEIF ( y < 0.0 )  THEN
92                alpha = 3.0 * pi / 2.0
93             ENDIF
94          ELSE
95             IF ( x < 0.0 )  THEN
96                alpha = ATAN( y / x ) + pi
97             ELSE
98                IF ( y < 0.0 )  THEN
99                   alpha = ATAN( y / x ) + 2.0 * pi
100                ELSE
101                   alpha = ATAN( y / x )
102                ENDIF
103             ENDIF
104          ENDIF
105
106          uw = -SIN( alpha ) * betrag
107
108          DO  k = kc1, kc2
109             u(k,j,i) = u(k,j,i) + uw
110          ENDDO
111       ENDDO
112    ENDDO
113
114!
115!-- Compute the v-component.
116    DO  i = nxl, nxr
117       DO  j = nys, nyn
118          x = ( i - ic ) * dx
119          y = ( j - jc - 0.5) * dy
120          radius = SQRT( x**2 + y**2 )
121          IF ( radius <= 2.0 * rc )  THEN
122             betrag = radius / ( 2.0 * rc ) * 0.08
123          ELSEIF ( radius > 2.0 * rc  .AND.  radius < 8.0 * rc )  THEN
124             betrag = 0.08 * EXP( -( radius - 2.0 * rc ) / 2.0 )
125          ELSE
126             betrag = 0.0
127          ENDIF
128
129          IF ( x == 0.0 )  THEN
130             IF ( y > 0.0 )  THEN
131                alpha = pi / 2.0
132             ELSEIF ( y < 0.0 )  THEN
133                alpha = 3.0 * pi / 2.0
134             ENDIF
135          ELSE
136             IF ( x < 0.0 )  THEN
137                alpha = ATAN( y / x ) + pi
138             ELSE
139                IF ( y < 0.0 )  THEN
140                   alpha = ATAN( y / x ) + 2.0 * pi
141                ELSE
142                   alpha = ATAN( y / x )
143                ENDIF
144             ENDIF
145          ENDIF
146
147          vw = COS( alpha ) * betrag
148
149          DO  k = kc1, kc2
150             v(k,j,i) = v(k,j,i) + vw
151          ENDDO
152       ENDDO
153    ENDDO
154
155!
156!-- Exchange of boundary values for the velocities.
157    CALL exchange_horiz( u, uxrp,    0 )
158    CALL exchange_horiz( v,    0, vynp )
159!
160!-- Make velocity field nondivergent.
161    n_sor = nsor_ini
162    CALL pres
163    n_sor = nsor
164
165 END SUBROUTINE init_rankine
Note: See TracBrowser for help on using the repository browser.