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

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

last commit documented

  • Property svn:keywords set to Id
File size: 5.0 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 1037 2012-10-22 14:10:22Z raasch $
26!
27! 1036 2012-10-22 13:43:42Z raasch
28! code put under GPL (PALM 3.9)
29!
30! 667 2010-12-23 12:06:00Z suehring/gryschka
31! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
32! Calls of exchange_horiz are modified.
33!
34! 107 2007-08-17 13:54:45Z raasch
35! Initial profiles are reset to constant profiles
36!
37! 75 2007-03-22 09:54:05Z raasch
38! uxrp, vynp eliminated, 2nd+3rd argument removed from exchange horiz
39!
40! RCS Log replace by Id keyword, revision history cleaned up
41!
42! Revision 1.11  2005/03/26 20:38:49  raasch
43! Arguments for non-cyclic boundary conditions added to argument list of
44! routine exchange_horiz
45!
46! Revision 1.1  1997/08/11 06:18:43  raasch
47! Initial revision
48!
49!
50! Description:
51! ------------
52! Initialize a (nondivergent) Rankine eddy with a vertical axis in order to test
53! the advection terms and the pressure solver.
54!------------------------------------------------------------------------------!
55
56    USE arrays_3d
57    USE constants
58    USE grid_variables
59    USE indices
60    USE control_parameters
61
62    IMPLICIT NONE
63
64    INTEGER ::  i, ic, j, jc, k, kc1, kc2
65    REAL    ::  alpha, betrag, radius, rc, uw, vw, x, y
66
67!
68!-- Default: eddy radius rc, eddy strength z,
69!--          position of eddy centre: ic, jc, kc1, kc2
70    rc  =  4.0 * dx
71    ic  =  ( nx+1 ) / 2
72    jc  =  ic
73    kc1 = nzb
74    kc2 = nzt+1
75
76!
77!-- Reset initial profiles to constant profiles
78    IF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 )  THEN
79       DO  i = nxlg, nxrg
80          DO  j = nysg, nyng
81             pt(:,j,i) = pt_init
82             u(:,j,i)  = u_init
83             v(:,j,i)  = v_init
84          ENDDO
85       ENDDO
86    ENDIF
87
88!
89!-- Compute the u-component.
90    DO  i = nxl, nxr
91       DO  j = nys, nyn
92          x = ( i - ic - 0.5 ) * dx
93          y = ( j - jc ) * 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          uw = -SIN( alpha ) * betrag
122
123          DO  k = kc1, kc2
124             u(k,j,i) = u(k,j,i) + uw
125          ENDDO
126       ENDDO
127    ENDDO
128
129!
130!-- Compute the v-component.
131    DO  i = nxl, nxr
132       DO  j = nys, nyn
133          x = ( i - ic ) * dx
134          y = ( j - jc - 0.5) * dy
135          radius = SQRT( x**2 + y**2 )
136          IF ( radius <= 2.0 * rc )  THEN
137             betrag = radius / ( 2.0 * rc ) * 0.08
138          ELSEIF ( radius > 2.0 * rc  .AND.  radius < 8.0 * rc )  THEN
139             betrag = 0.08 * EXP( -( radius - 2.0 * rc ) / 2.0 )
140          ELSE
141             betrag = 0.0
142          ENDIF
143
144          IF ( x == 0.0 )  THEN
145             IF ( y > 0.0 )  THEN
146                alpha = pi / 2.0
147             ELSEIF ( y < 0.0 )  THEN
148                alpha = 3.0 * pi / 2.0
149             ENDIF
150          ELSE
151             IF ( x < 0.0 )  THEN
152                alpha = ATAN( y / x ) + pi
153             ELSE
154                IF ( y < 0.0 )  THEN
155                   alpha = ATAN( y / x ) + 2.0 * pi
156                ELSE
157                   alpha = ATAN( y / x )
158                ENDIF
159             ENDIF
160          ENDIF
161
162          vw = COS( alpha ) * betrag
163
164          DO  k = kc1, kc2
165             v(k,j,i) = v(k,j,i) + vw
166          ENDDO
167       ENDDO
168    ENDDO
169
170!
171!-- Exchange of boundary values for the velocities.
172    CALL exchange_horiz( u, nbgp)
173    CALL exchange_horiz( v, nbgp )
174!
175!-- Make velocity field nondivergent.
176    n_sor = nsor_ini
177    CALL pres
178    n_sor = nsor
179
180 END SUBROUTINE init_rankine
Note: See TracBrowser for help on using the repository browser.