source: palm/trunk/SOURCE/coriolis.f90 @ 4679

Last change on this file since 4679 was 4559, checked in by raasch, 4 years ago

files re-formatted to follow the PALM coding standard

  • Property svn:keywords set to Id
File size: 12.0 KB
Line 
1!> @file coriolis.f90
2!--------------------------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
6! Public License as published by the Free Software Foundation, either version 3 of the License, or
7! (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11! Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with PALM. If not, see
14! <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2020 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: coriolis.f90 4559 2020-06-11 08:51:48Z pavelkrc $
26! file re-formatted to follow the PALM coding standard
27!
28! 4360 2020-01-07 11:25:50Z suehring
29! Introduction of wall_flags_total_0, which currently sets bits based on static topography
30! information used in wall_flags_static_0
31!
32! 4329 2019-12-10 15:46:36Z motisi
33! Renamed wall_flags_0 to wall_flags_static_0
34!
35! 4196 2019-08-29 11:02:06Z gronemeier
36! Consider rotation of model domain
37!
38! 4182 2019-08-22 15:20:23Z scharf
39! Corrected "Former revisions" section
40!
41! 3655 2019-01-07 16:51:22Z knoop
42! OpenACC port for SPEC
43!
44! Revision 1.1  1997/08/29 08:57:38  raasch
45! Initial revision
46!
47!
48! Description:
49! ------------
50!> Computation of all Coriolis terms in the equations of motion.
51!>
52!> @note In this routine the topography is masked, even though this is again done in
53!> prognostic_equations. However, omitting the masking here lead to slightly different results.
54!> Reason unknown.
55!--------------------------------------------------------------------------------------------------!
56 MODULE coriolis_mod
57
58
59    PRIVATE
60    PUBLIC coriolis
61
62    INTERFACE coriolis
63       MODULE PROCEDURE coriolis
64       MODULE PROCEDURE coriolis_ij
65    END INTERFACE coriolis
66
67 CONTAINS
68
69
70!--------------------------------------------------------------------------------------------------!
71! Description:
72! ------------
73!> Call for all grid points
74!--------------------------------------------------------------------------------------------------!
75    SUBROUTINE coriolis( component )
76
77       USE arrays_3d,                                                                              &
78           ONLY:  tend, u, ug, v, vg, w
79
80       USE basic_constants_and_equations_mod,                                                      &
81           ONLY:  pi
82
83       USE control_parameters,                                                                     &
84           ONLY:  f, fs, message_string, rotation_angle
85
86       USE indices,                                                                                &
87           ONLY:  nxl, nxlu, nxr, nyn, nys, nysv, nzb, nzt, wall_flags_total_0
88
89       USE kinds
90
91       IMPLICIT NONE
92
93       INTEGER(iwp) ::  component      !< component of momentum equation
94       INTEGER(iwp) ::  i              !< running index x direction
95       INTEGER(iwp) ::  j              !< running index y direction
96       INTEGER(iwp) ::  k              !< running index z direction
97
98       REAL(wp)     ::  cos_rot_angle  !< cosine of model rotation angle
99       REAL(wp)     ::  flag           !< flag to mask topography
100       REAL(wp)     ::  sin_rot_angle  !< sine of model rotation angle
101
102!
103!--    Precalculate cosine and sine of rotation angle
104       cos_rot_angle = COS( rotation_angle * pi / 180.0_wp )
105       sin_rot_angle = SIN( rotation_angle * pi / 180.0_wp )
106
107!
108!--    Compute Coriolis terms for the three velocity components
109       SELECT CASE ( component )
110
111!
112!--       u-component
113          CASE ( 1 )
114             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k, flag) &
115             !$ACC PRESENT(wall_flags_total_0) &
116             !$ACC PRESENT(v, w, vg) &
117             !$ACC PRESENT(tend)
118             DO  i = nxlu, nxr
119                DO  j = nys, nyn
120                   DO  k = nzb+1, nzt
121!
122!--                   Predetermine flag to mask topography
123                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
124
125                      tend(k,j,i) = tend(k,j,i) + flag *                                           &
126                            ( f                                                                    &
127                              * ( 0.25_wp * ( v(k,j,i-1) + v(k,j,i) + v(k,j+1,i-1) + v(k,j+1,i) )  &
128                                - vg(k) )                                                          &
129                            - fs * cos_rot_angle                                                   &
130                              * 0.25_wp * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) )    &
131                            )
132                   ENDDO
133                ENDDO
134             ENDDO
135
136!
137!--       v-component
138          CASE ( 2 )
139             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k, flag) &
140             !$ACC PRESENT(wall_flags_total_0) &
141             !$ACC PRESENT(u, w, ug) &
142             !$ACC PRESENT(tend)
143             DO  i = nxl, nxr
144                DO  j = nysv, nyn
145                   DO  k = nzb+1, nzt
146!
147!--                   Predetermine flag to mask topography
148                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
149
150                      tend(k,j,i) = tend(k,j,i) - flag *                                           &
151                            ( f                                                                    &
152                              * ( 0.25_wp * ( u(k,j-1,i) + u(k,j,i) + u(k,j-1,i+1) + u(k,j,i+1) )  &
153                                - ug(k) )                                                          &
154                            + fs * sin_rot_angle                                                   &
155                              * 0.25_wp * ( w(k,j,i) + w(k-1,j,i) + w(k,j-1,i) + w(k-1,j-1,i) )    &
156                            )
157                   ENDDO
158                ENDDO
159             ENDDO
160
161!
162!--       w-component
163          CASE ( 3 )
164             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k, flag) &
165             !$ACC PRESENT(wall_flags_total_0) &
166             !$ACC PRESENT(u, v) &
167             !$ACC PRESENT(tend)
168             DO  i = nxl, nxr
169                DO  j = nys, nyn
170                   DO  k = nzb+1, nzt
171!
172!--                   Predetermine flag to mask topography
173                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
174
175                      tend(k,j,i) = tend(k,j,i)                                                    &
176                                     + fs * 0.25_wp * flag                                         &
177                                       * ( cos_rot_angle                                           &
178                                           * ( u(k,j,i) + u(k+1,j,i) + u(k,j,i+1) + u(k+1,j,i+1) ) &
179                                         + sin_rot_angle                                           &
180                                           * ( v(k,j,i) + v(k+1,j,i) + v(k,j+1,i) + v(k+1,j+1,i) ) &
181                                         )
182                   ENDDO
183                ENDDO
184             ENDDO
185
186          CASE DEFAULT
187
188             WRITE( message_string, * ) ' wrong component: ', component
189             CALL message( 'coriolis', 'PA0173', 1, 2, 0, 6, 0 )
190
191       END SELECT
192
193    END SUBROUTINE coriolis
194
195
196!--------------------------------------------------------------------------------------------------!
197! Description:
198! ------------
199!> Call for grid point i,j
200!--------------------------------------------------------------------------------------------------!
201    SUBROUTINE coriolis_ij( i, j, component )
202
203       USE arrays_3d,                                                                              &
204           ONLY:  tend, u, ug, v, vg, w
205
206       USE basic_constants_and_equations_mod,                                                      &
207           ONLY:  pi
208
209       USE control_parameters,                                                                     &
210           ONLY:  f, fs, message_string, rotation_angle
211
212       USE indices,                                                                                &
213           ONLY:  nzb, nzt, wall_flags_total_0
214
215       USE kinds
216
217       IMPLICIT NONE
218
219       INTEGER(iwp) ::  component  !< component of momentum equation
220       INTEGER(iwp) ::  i          !< running index x direction
221       INTEGER(iwp) ::  j          !< running index y direction
222       INTEGER(iwp) ::  k          !< running index z direction
223
224       REAL(wp)     ::  cos_rot_angle  !< cosine of model rotation angle
225       REAL(wp)     ::  flag           !< flag to mask topography
226       REAL(wp)     ::  sin_rot_angle  !< sine of model rotation angle
227
228!
229!--    Precalculate cosine and sine of rotation angle
230       cos_rot_angle = COS( rotation_angle * pi / 180.0_wp )
231       sin_rot_angle = SIN( rotation_angle * pi / 180.0_wp )
232
233!
234!--    Compute Coriolis terms for the three velocity components
235       SELECT CASE ( component )
236
237!
238!--       u-component
239          CASE ( 1 )
240             DO  k = nzb+1, nzt
241!
242!--             Predetermine flag to mask topography
243                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
244
245                tend(k,j,i) = tend(k,j,i) + flag *                                                 &
246                            ( f                                                                    &
247                              * ( 0.25_wp * ( v(k,j,i-1) + v(k,j,i) + v(k,j+1,i-1) + v(k,j+1,i) )  &
248                                - vg(k) )                                                          &
249                            - fs * cos_rot_angle                                                   &
250                              * 0.25_wp * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) )    &
251                            )
252             ENDDO
253
254!
255!--       v-component
256          CASE ( 2 )
257             DO  k = nzb+1, nzt
258!
259!--             Predetermine flag to mask topography
260                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
261
262                tend(k,j,i) = tend(k,j,i) - flag *                                                 &
263                              ( f                                                                  &
264                                * ( 0.25_wp * ( u(k,j-1,i) + u(k,j,i) + u(k,j-1,i+1) + u(k,j,i+1) )&
265                                  - ug(k) )                                                        &
266                              + fs * sin_rot_angle                                                 &
267                                * 0.25_wp * ( w(k,j,i) + w(k-1,j,i) + w(k,j-1,i) + w(k-1,j-1,i) )  &
268                              )
269             ENDDO
270
271!
272!--       w-component
273          CASE ( 3 )
274             DO  k = nzb+1, nzt
275!
276!--             Predetermine flag to mask topography
277                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
278
279                tend(k,j,i) = tend(k,j,i)                                                          &
280                              + fs * 0.25_wp * flag                                                &
281                                * ( cos_rot_angle                                                  &
282                                    * ( u(k,j,i) + u(k+1,j,i) + u(k,j,i+1) + u(k+1,j,i+1) )        &
283                                  + sin_rot_angle                                                  &
284                                    * ( v(k,j,i) + v(k+1,j,i) + v(k,j+1,i) + v(k+1,j+1,i) )        &
285                                  )
286             ENDDO
287
288          CASE DEFAULT
289
290             WRITE( message_string, * ) ' wrong component: ', component
291             CALL message( 'coriolis', 'PA0173', 1, 2, 0, 6, 0 )
292
293       END SELECT
294
295    END SUBROUTINE coriolis_ij
296
297 END MODULE coriolis_mod
Note: See TracBrowser for help on using the repository browser.