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

Last change on this file since 1354 was 1354, checked in by heinze, 7 years ago

last commit documented

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