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

Last change on this file since 1322 was 1322, checked in by raasch, 8 years ago

REAL functions and a lot of REAL constants provided with KIND-attribute,
some small bugfixes

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