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

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