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

Last change on this file since 1036 was 1036, checked in by raasch, 9 years ago

code has been put under the GNU General Public License (v3)

  • Property svn:keywords set to Id
File size: 4.9 KB
Line 
1 SUBROUTINE init_rankine
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2012  Leibniz University Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23! Former revisions:
24! -----------------
25! $Id: init_rankine.f90 1036 2012-10-22 13:43:42Z raasch $
26!
27! 667 2010-12-23 12:06:00Z suehring/gryschka
28! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
29! Calls of exchange_horiz are modified.
30!
31! 107 2007-08-17 13:54:45Z raasch
32! Initial profiles are reset to constant profiles
33!
34! 75 2007-03-22 09:54:05Z raasch
35! uxrp, vynp eliminated, 2nd+3rd argument removed from exchange horiz
36!
37! RCS Log replace by Id keyword, revision history cleaned up
38!
39! Revision 1.11  2005/03/26 20:38:49  raasch
40! Arguments for non-cyclic boundary conditions added to argument list of
41! routine exchange_horiz
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!-- Reset initial profiles to constant profiles
75    IF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 )  THEN
76       DO  i = nxlg, nxrg
77          DO  j = nysg, nyng
78             pt(:,j,i) = pt_init
79             u(:,j,i)  = u_init
80             v(:,j,i)  = v_init
81          ENDDO
82       ENDDO
83    ENDIF
84
85!
86!-- Compute the u-component.
87    DO  i = nxl, nxr
88       DO  j = nys, nyn
89          x = ( i - ic - 0.5 ) * dx
90          y = ( j - jc ) * dy
91          radius = SQRT( x**2 + y**2 )
92          IF ( radius <= 2.0 * rc )  THEN
93             betrag = radius / ( 2.0 * rc ) * 0.08
94          ELSEIF ( radius > 2.0 * rc  .AND.  radius < 8.0 * rc )  THEN
95             betrag = 0.08 * EXP( -( radius - 2.0 * rc ) / 2.0 )
96          ELSE
97             betrag = 0.0
98          ENDIF
99
100          IF ( x == 0.0 )  THEN
101             IF ( y > 0.0 )  THEN
102                alpha = pi / 2.0
103             ELSEIF ( y < 0.0 )  THEN
104                alpha = 3.0 * pi / 2.0
105             ENDIF
106          ELSE
107             IF ( x < 0.0 )  THEN
108                alpha = ATAN( y / x ) + pi
109             ELSE
110                IF ( y < 0.0 )  THEN
111                   alpha = ATAN( y / x ) + 2.0 * pi
112                ELSE
113                   alpha = ATAN( y / x )
114                ENDIF
115             ENDIF
116          ENDIF
117
118          uw = -SIN( alpha ) * betrag
119
120          DO  k = kc1, kc2
121             u(k,j,i) = u(k,j,i) + uw
122          ENDDO
123       ENDDO
124    ENDDO
125
126!
127!-- Compute the v-component.
128    DO  i = nxl, nxr
129       DO  j = nys, nyn
130          x = ( i - ic ) * dx
131          y = ( j - jc - 0.5) * dy
132          radius = SQRT( x**2 + y**2 )
133          IF ( radius <= 2.0 * rc )  THEN
134             betrag = radius / ( 2.0 * rc ) * 0.08
135          ELSEIF ( radius > 2.0 * rc  .AND.  radius < 8.0 * rc )  THEN
136             betrag = 0.08 * EXP( -( radius - 2.0 * rc ) / 2.0 )
137          ELSE
138             betrag = 0.0
139          ENDIF
140
141          IF ( x == 0.0 )  THEN
142             IF ( y > 0.0 )  THEN
143                alpha = pi / 2.0
144             ELSEIF ( y < 0.0 )  THEN
145                alpha = 3.0 * pi / 2.0
146             ENDIF
147          ELSE
148             IF ( x < 0.0 )  THEN
149                alpha = ATAN( y / x ) + pi
150             ELSE
151                IF ( y < 0.0 )  THEN
152                   alpha = ATAN( y / x ) + 2.0 * pi
153                ELSE
154                   alpha = ATAN( y / x )
155                ENDIF
156             ENDIF
157          ENDIF
158
159          vw = COS( alpha ) * betrag
160
161          DO  k = kc1, kc2
162             v(k,j,i) = v(k,j,i) + vw
163          ENDDO
164       ENDDO
165    ENDDO
166
167!
168!-- Exchange of boundary values for the velocities.
169    CALL exchange_horiz( u, nbgp)
170    CALL exchange_horiz( v, nbgp )
171!
172!-- Make velocity field nondivergent.
173    n_sor = nsor_ini
174    CALL pres
175    n_sor = nsor
176
177 END SUBROUTINE init_rankine
Note: See TracBrowser for help on using the repository browser.