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

Last change on this file since 4765 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
RevLine 
[1873]1!> @file coriolis.f90
[4559]2!--------------------------------------------------------------------------------------------------!
[2696]3! This file is part of the PALM model system.
[1036]4!
[4559]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.
[1036]8!
[4559]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.
[1036]12!
[4559]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/>.
[1036]15!
[4360]16! Copyright 1997-2020 Leibniz Universitaet Hannover
[4559]17!--------------------------------------------------------------------------------------------------!
[1036]18!
[254]19! Current revisions:
[1]20! -----------------
[1354]21!
[3183]22!
[1321]23! Former revisions:
24! -----------------
25! $Id: coriolis.f90 4559 2020-06-11 08:51:48Z Giersch $
[4559]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!
[4346]32! 4329 2019-12-10 15:46:36Z motisi
[4329]33! Renamed wall_flags_0 to wall_flags_static_0
[4559]34!
[4329]35! 4196 2019-08-29 11:02:06Z gronemeier
[4196]36! Consider rotation of model domain
[4559]37!
[4196]38! 4182 2019-08-22 15:20:23Z scharf
[4182]39! Corrected "Former revisions" section
[4559]40!
[4182]41! 3655 2019-01-07 16:51:22Z knoop
[3634]42! OpenACC port for SPEC
[1321]43!
[4182]44! Revision 1.1  1997/08/29 08:57:38  raasch
45! Initial revision
46!
47!
[1]48! Description:
49! ------------
[1682]50!> Computation of all Coriolis terms in the equations of motion.
[4559]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!--------------------------------------------------------------------------------------------------!
[1682]56 MODULE coriolis_mod
[1]57
[4559]58
[1]59    PRIVATE
[2118]60    PUBLIC coriolis
[1]61
62    INTERFACE coriolis
63       MODULE PROCEDURE coriolis
64       MODULE PROCEDURE coriolis_ij
65    END INTERFACE coriolis
66
67 CONTAINS
68
69
[4559]70!--------------------------------------------------------------------------------------------------!
[1682]71! Description:
72! ------------
73!> Call for all grid points
[4559]74!--------------------------------------------------------------------------------------------------!
[1]75    SUBROUTINE coriolis( component )
76
[4559]77       USE arrays_3d,                                                                              &
78           ONLY:  tend, u, ug, v, vg, w
[4196]79
[4559]80       USE basic_constants_and_equations_mod,                                                      &
[4196]81           ONLY:  pi
82
[4559]83       USE control_parameters,                                                                     &
[4196]84           ONLY:  f, fs, message_string, rotation_angle
[4559]85
86       USE indices,                                                                                &
87           ONLY:  nxl, nxlu, nxr, nyn, nys, nysv, nzb, nzt, wall_flags_total_0
88
[1320]89       USE kinds
[1]90
91       IMPLICIT NONE
92
[4196]93       INTEGER(iwp) ::  component      !< component of momentum equation
[4559]94       INTEGER(iwp) ::  i              !< running index x direction
95       INTEGER(iwp) ::  j              !< running index y direction
96       INTEGER(iwp) ::  k              !< running index z direction
[1]97
[4196]98       REAL(wp)     ::  cos_rot_angle  !< cosine of model rotation angle
[3182]99       REAL(wp)     ::  flag           !< flag to mask topography
[4196]100       REAL(wp)     ::  sin_rot_angle  !< sine of model rotation angle
[2696]101
[1]102!
[4196]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!
[1]108!--    Compute Coriolis terms for the three velocity components
109       SELECT CASE ( component )
110
111!
112!--       u-component
113          CASE ( 1 )
[3634]114             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k, flag) &
[4346]115             !$ACC PRESENT(wall_flags_total_0) &
[3634]116             !$ACC PRESENT(v, w, vg) &
117             !$ACC PRESENT(tend)
[106]118             DO  i = nxlu, nxr
[1]119                DO  j = nys, nyn
[2232]120                   DO  k = nzb+1, nzt
121!
[2696]122!--                   Predetermine flag to mask topography
[4346]123                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
[2232]124
[4196]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                            )
[1]132                   ENDDO
133                ENDDO
134             ENDDO
135
136!
137!--       v-component
138          CASE ( 2 )
[3634]139             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k, flag) &
[4346]140             !$ACC PRESENT(wall_flags_total_0) &
[4196]141             !$ACC PRESENT(u, w, ug) &
[3634]142             !$ACC PRESENT(tend)
[1]143             DO  i = nxl, nxr
[106]144                DO  j = nysv, nyn
[2232]145                   DO  k = nzb+1, nzt
[2696]146!
147!--                   Predetermine flag to mask topography
[4346]148                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
[2696]149
[4196]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                            )
[1]157                   ENDDO
158                ENDDO
159             ENDDO
160
161!
162!--       w-component
163          CASE ( 3 )
[3634]164             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k, flag) &
[4346]165             !$ACC PRESENT(wall_flags_total_0) &
[4196]166             !$ACC PRESENT(u, v) &
[3634]167             !$ACC PRESENT(tend)
[1]168             DO  i = nxl, nxr
169                DO  j = nys, nyn
[2232]170                   DO  k = nzb+1, nzt
[2696]171!
172!--                   Predetermine flag to mask topography
[4346]173                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
[2696]174
[4559]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                                         )
[1]182                   ENDDO
183                ENDDO
184             ENDDO
185
186          CASE DEFAULT
187
[254]188             WRITE( message_string, * ) ' wrong component: ', component
189             CALL message( 'coriolis', 'PA0173', 1, 2, 0, 6, 0 )
[1]190
191       END SELECT
192
193    END SUBROUTINE coriolis
194
195
[4559]196!--------------------------------------------------------------------------------------------------!
[1682]197! Description:
198! ------------
199!> Call for grid point i,j
[4559]200!--------------------------------------------------------------------------------------------------!
[1]201    SUBROUTINE coriolis_ij( i, j, component )
202
[4559]203       USE arrays_3d,                                                                              &
204           ONLY:  tend, u, ug, v, vg, w
[4196]205
[4559]206       USE basic_constants_and_equations_mod,                                                      &
[4196]207           ONLY:  pi
208
[4559]209       USE control_parameters,                                                                     &
[4196]210           ONLY:  f, fs, message_string, rotation_angle
[4559]211
212       USE indices,                                                                                &
[4346]213           ONLY:  nzb, nzt, wall_flags_total_0
[4559]214
[1320]215       USE kinds
[4196]216
[1]217       IMPLICIT NONE
218
[4196]219       INTEGER(iwp) ::  component  !< component of momentum equation
[4559]220       INTEGER(iwp) ::  i          !< running index x direction
221       INTEGER(iwp) ::  j          !< running index y direction
222       INTEGER(iwp) ::  k          !< running index z direction
[1]223
[4196]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
[2232]227
[1]228!
[4196]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!
[1]234!--    Compute Coriolis terms for the three velocity components
235       SELECT CASE ( component )
236
237!
238!--       u-component
239          CASE ( 1 )
[2232]240             DO  k = nzb+1, nzt
241!
242!--             Predetermine flag to mask topography
[4346]243                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
[2232]244
[4196]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                            )
[1]252             ENDDO
253
254!
255!--       v-component
256          CASE ( 2 )
[2232]257             DO  k = nzb+1, nzt
[2696]258!
259!--             Predetermine flag to mask topography
[4346]260                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
[2696]261
[4196]262                tend(k,j,i) = tend(k,j,i) - flag *                                                 &
[4559]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                              )
[1]269             ENDDO
270
271!
272!--       w-component
273          CASE ( 3 )
[2232]274             DO  k = nzb+1, nzt
[2696]275!
276!--             Predetermine flag to mask topography
[4346]277                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
[2696]278
[4559]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                                  )
[1]286             ENDDO
287
288          CASE DEFAULT
289
[254]290             WRITE( message_string, * ) ' wrong component: ', component
291             CALL message( 'coriolis', 'PA0173', 1, 2, 0, 6, 0 )
[1]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.