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

Last change on this file since 1321 was 1321, checked in by raasch, 7 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 5.7 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-2014 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23! Former revisions:
24! -----------------
25! $Id: init_rankine.f90 1321 2014-03-20 09:40:40Z raasch $
26!
27! 1320 2014-03-20 08:40:49Z raasch
28! ONLY-attribute added to USE-statements,
29! kind-parameters added to all INTEGER and REAL declaration statements,
30! kinds are defined in new module kinds,
31! revision history before 2012 removed,
32! comment fields (!:) to be used for variable explanations added to
33! all variable declaration statements
34!
35! 1036 2012-10-22 13:43:42Z raasch
36! code put under GPL (PALM 3.9)
37!
38! Revision 1.1  1997/08/11 06:18:43  raasch
39! Initial revision
40!
41!
42! Description:
43! ------------
44! Initialize a (nondivergent) Rankine eddy with a vertical axis in order to test
45! the advection terms and the pressure solver.
46!------------------------------------------------------------------------------!
47
48    USE arrays_3d,                                                             &
49        ONLY:  pt, pt_init, u, u_init, v, v_init
50
51    USE control_parameters,                                                    &
52        ONLY:  initializing_actions, n_sor, nsor, nsor_ini   
53
54    USE constants,                                                             &
55        ONLY:  pi
56
57    USE grid_variables,                                                        &
58        ONLY:  dx, dy
59
60    USE indices,                                                               &
61        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt     
62               
63    USE kinds
64
65    IMPLICIT NONE
66
67    INTEGER(iwp) ::  i   !:
68    INTEGER(iwp) ::  ic  !:
69    INTEGER(iwp) ::  j   !:
70    INTEGER(iwp) ::  jc  !:
71    INTEGER(iwp) ::  k   !:
72    INTEGER(iwp) ::  kc1 !:
73    INTEGER(iwp) ::  kc2 !:
74   
75    REAL(wp)     ::  alpha  !:
76    REAL(wp)     ::  betrag !:
77    REAL(wp)     ::  radius !:
78    REAL(wp)     ::  rc     !:
79    REAL(wp)     ::  uw     !:
80    REAL(wp)     ::  vw     !:
81    REAL(wp)     ::  x      !:
82    REAL(wp)     ::  y      !:
83
84!
85!-- Default: eddy radius rc, eddy strength z,
86!--          position of eddy centre: ic, jc, kc1, kc2
87    rc  =  4.0 * dx
88    ic  =  ( nx+1 ) / 2
89    jc  =  ic
90    kc1 = nzb
91    kc2 = nzt+1
92
93!
94!-- Reset initial profiles to constant profiles
95    IF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 )  THEN
96       DO  i = nxlg, nxrg
97          DO  j = nysg, nyng
98             pt(:,j,i) = pt_init
99             u(:,j,i)  = u_init
100             v(:,j,i)  = v_init
101          ENDDO
102       ENDDO
103    ENDIF
104
105!
106!-- Compute the u-component.
107    DO  i = nxl, nxr
108       DO  j = nys, nyn
109          x = ( i - ic - 0.5 ) * dx
110          y = ( j - jc ) * dy
111          radius = SQRT( x**2 + y**2 )
112          IF ( radius <= 2.0 * rc )  THEN
113             betrag = radius / ( 2.0 * rc ) * 0.08
114          ELSEIF ( radius > 2.0 * rc  .AND.  radius < 8.0 * rc )  THEN
115             betrag = 0.08 * EXP( -( radius - 2.0 * rc ) / 2.0 )
116          ELSE
117             betrag = 0.0
118          ENDIF
119
120          IF ( x == 0.0 )  THEN
121             IF ( y > 0.0 )  THEN
122                alpha = pi / 2.0
123             ELSEIF ( y < 0.0 )  THEN
124                alpha = 3.0 * pi / 2.0
125             ENDIF
126          ELSE
127             IF ( x < 0.0 )  THEN
128                alpha = ATAN( y / x ) + pi
129             ELSE
130                IF ( y < 0.0 )  THEN
131                   alpha = ATAN( y / x ) + 2.0 * pi
132                ELSE
133                   alpha = ATAN( y / x )
134                ENDIF
135             ENDIF
136          ENDIF
137
138          uw = -SIN( alpha ) * betrag
139
140          DO  k = kc1, kc2
141             u(k,j,i) = u(k,j,i) + uw
142          ENDDO
143       ENDDO
144    ENDDO
145
146!
147!-- Compute the v-component.
148    DO  i = nxl, nxr
149       DO  j = nys, nyn
150          x = ( i - ic ) * dx
151          y = ( j - jc - 0.5) * dy
152          radius = SQRT( x**2 + y**2 )
153          IF ( radius <= 2.0 * rc )  THEN
154             betrag = radius / ( 2.0 * rc ) * 0.08
155          ELSEIF ( radius > 2.0 * rc  .AND.  radius < 8.0 * rc )  THEN
156             betrag = 0.08 * EXP( -( radius - 2.0 * rc ) / 2.0 )
157          ELSE
158             betrag = 0.0
159          ENDIF
160
161          IF ( x == 0.0 )  THEN
162             IF ( y > 0.0 )  THEN
163                alpha = pi / 2.0
164             ELSEIF ( y < 0.0 )  THEN
165                alpha = 3.0 * pi / 2.0
166             ENDIF
167          ELSE
168             IF ( x < 0.0 )  THEN
169                alpha = ATAN( y / x ) + pi
170             ELSE
171                IF ( y < 0.0 )  THEN
172                   alpha = ATAN( y / x ) + 2.0 * pi
173                ELSE
174                   alpha = ATAN( y / x )
175                ENDIF
176             ENDIF
177          ENDIF
178
179          vw = COS( alpha ) * betrag
180
181          DO  k = kc1, kc2
182             v(k,j,i) = v(k,j,i) + vw
183          ENDDO
184       ENDDO
185    ENDDO
186
187!
188!-- Exchange of boundary values for the velocities.
189    CALL exchange_horiz( u, nbgp)
190    CALL exchange_horiz( v, nbgp )
191!
192!-- Make velocity field nondivergent.
193    n_sor = nsor_ini
194    CALL pres
195    n_sor = nsor
196
197 END SUBROUTINE init_rankine
Note: See TracBrowser for help on using the repository browser.