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

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

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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