Home

Context Navigation

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

Last change on this file since 4304 was 4196, checked in by gronemeier, 5 years ago

Consider rotation of model domain for claculation of Coriolis force:

check_parameters: Overwrite rotation_angle from namelist by value from static driver;
coriolis: Consider rotation of model domain;
header: Write information about rotation angle;
modules: Added rotation_angle;
netcdf_interface_mod: replaced rotation angle from input-netCDF file by namelist parameter 'rotation_angle';
ocean_mod: Consider rotation of model domain for calculating the Stokes drift;
parin: added rotation_angle to initialization_parameters;

Property svn:keywords set to Id

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |