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

Last change on this file since 1032 was 668, checked in by suehring, 14 years ago

last commit documented

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