Home

Context Navigation

source: palm/trunk/SOURCE/diffusion_s.f90 @ 3634

Last change on this file since 3634 was 3634, checked in by knoop, 5 years ago
OpenACC port for SPEC
Property svn:keywords set to `Id`
File size: 34.7 KB

Line
1	!> @file diffusion_s.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-2018 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: diffusion_s.f90 3634 2018-12-18 12:31:28Z knoop $
27	! OpenACC port for SPEC
28	!
29	! 3547 2018-11-21 13:21:24Z suehring
30	! variables documented
31	!
32	! 2759 2018-01-17 16:24:59Z suehring
33	! Major bugfix, horizontal diffusion at vertical surfaces corrected.
34	!
35	! 2718 2018-01-02 08:49:38Z maronga
36	! Corrected "Former revisions" section
37	!
38	! 2696 2017-12-14 17:12:51Z kanani
39	! Change in file header (GPL part)
40	!
41	! 2233 2017-05-30 18:08:54Z suehring
42	!
43	! 2232 2017-05-30 17:47:52Z suehring
44	! Adjustments to new topography and surface concept
45	!
46	! 2118 2017-01-17 16:38:49Z raasch
47	! OpenACC version of subroutine removed
48	!
49	! 2037 2016-10-26 11:15:40Z knoop
50	! Anelastic approximation implemented
51	!
52	! 2000 2016-08-20 18:09:15Z knoop
53	! Forced header and separation lines into 80 columns
54	!
55	! 1873 2016-04-18 14:50:06Z maronga
56	! Module renamed (removed _mod)
57	!
58	! 1850 2016-04-08 13:29:27Z maronga
59	! Module renamed
60	!
61	! 1691 2015-10-26 16:17:44Z maronga
62	! Formatting corrections.
63	!
64	! 1682 2015-10-07 23:56:08Z knoop
65	! Code annotations made doxygen readable
66	!
67	! 1374 2014-04-25 12:55:07Z raasch
68	! missing variables added to ONLY list
69	!
70	! 1340 2014-03-25 19:45:13Z kanani
71	! REAL constants defined as wp-kind
72	!
73	! 1320 2014-03-20 08:40:49Z raasch
74	! ONLY-attribute added to USE-statements,
75	! kind-parameters added to all INTEGER and REAL declaration statements,
76	! kinds are defined in new module kinds,
77	! revision history before 2012 removed,
78	! comment fields (!:) to be used for variable explanations added to
79	! all variable declaration statements
80	!
81	! 1257 2013-11-08 15:18:40Z raasch
82	! openacc loop and loop vector clauses removed
83	!
84	! 1128 2013-04-12 06:19:32Z raasch
85	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
86	! j_north
87	!
88	! 1092 2013-02-02 11:24:22Z raasch
89	! unused variables removed
90	!
91	! 1036 2012-10-22 13:43:42Z raasch
92	! code put under GPL (PALM 3.9)
93	!
94	! 1015 2012-09-27 09:23:24Z raasch
95	! accelerator version (*_acc) added
96	!
97	! 1010 2012-09-20 07:59:54Z raasch
98	! cpp switch __nopointer added for pointer free version
99	!
100	! 1001 2012-09-13 14:08:46Z raasch
101	! some arrays comunicated by module instead of parameter list
102	!
103	! Revision 1.1 2000/04/13 14:54:02 schroeter
104	! Initial revision
105	!
106	!
107	! Description:
108	! ------------
109	!> Diffusion term of scalar quantities (temperature and water content)
110	!------------------------------------------------------------------------------!
111	MODULE diffusion_s_mod
112
113
114	PRIVATE
115	PUBLIC diffusion_s
116
117	INTERFACE diffusion_s
118	MODULE PROCEDURE diffusion_s
119	MODULE PROCEDURE diffusion_s_ij
120	END INTERFACE diffusion_s
121
122	CONTAINS
123
124
125	!------------------------------------------------------------------------------!
126	! Description:
127	! ------------
128	!> Call for all grid points
129	!------------------------------------------------------------------------------!
130	SUBROUTINE diffusion_s( s, s_flux_def_h_up, s_flux_def_h_down, &
131	s_flux_t, &
132	s_flux_lsm_h_up, s_flux_usm_h_up, &
133	s_flux_def_v_north, s_flux_def_v_south, &
134	s_flux_def_v_east, s_flux_def_v_west, &
135	s_flux_lsm_v_north, s_flux_lsm_v_south, &
136	s_flux_lsm_v_east, s_flux_lsm_v_west, &
137	s_flux_usm_v_north, s_flux_usm_v_south, &
138	s_flux_usm_v_east, s_flux_usm_v_west )
139
140	USE arrays_3d, &
141	ONLY: ddzu, ddzw, kh, tend, drho_air, rho_air_zw
142
143	USE control_parameters, &
144	ONLY: use_surface_fluxes, use_top_fluxes
145
146	USE grid_variables, &
147	ONLY: ddx, ddx2, ddy, ddy2
148
149	USE indices, &
150	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, &
151	nzt, wall_flags_0
152
153	USE kinds
154
155	USE surface_mod, &
156	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
157	surf_usm_v
158
159	IMPLICIT NONE
160
161	INTEGER(iwp) :: i !< running index x direction
162	INTEGER(iwp) :: j !< running index y direction
163	INTEGER(iwp) :: k !< running index z direction
164	INTEGER(iwp) :: m !< running index surface elements
165	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
166	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
167
168	REAL(wp) :: flag !< flag to mask topography grid points
169	REAL(wp) :: mask_bottom !< flag to mask vertical upward-facing surface
170	REAL(wp) :: mask_east !< flag to mask vertical surface east of the grid point
171	REAL(wp) :: mask_north !< flag to mask vertical surface north of the grid point
172	REAL(wp) :: mask_south !< flag to mask vertical surface south of the grid point
173	REAL(wp) :: mask_west !< flag to mask vertical surface west of the grid point
174	REAL(wp) :: mask_top !< flag to mask vertical downward-facing surface
175
176	REAL(wp), DIMENSION(1:surf_def_v(0)%ns) :: s_flux_def_v_north !< flux at north-facing vertical default-type surfaces
177	REAL(wp), DIMENSION(1:surf_def_v(1)%ns) :: s_flux_def_v_south !< flux at south-facing vertical default-type surfaces
178	REAL(wp), DIMENSION(1:surf_def_v(2)%ns) :: s_flux_def_v_east !< flux at east-facing vertical default-type surfaces
179	REAL(wp), DIMENSION(1:surf_def_v(3)%ns) :: s_flux_def_v_west !< flux at west-facing vertical default-type surfaces
180	REAL(wp), DIMENSION(1:surf_def_h(0)%ns) :: s_flux_def_h_up !< flux at horizontal upward-facing default-type surfaces
181	REAL(wp), DIMENSION(1:surf_def_h(1)%ns) :: s_flux_def_h_down !< flux at horizontal donwward-facing default-type surfaces
182	REAL(wp), DIMENSION(1:surf_lsm_h%ns) :: s_flux_lsm_h_up !< flux at horizontal upward-facing natural-type surfaces
183	REAL(wp), DIMENSION(1:surf_lsm_v(0)%ns) :: s_flux_lsm_v_north !< flux at north-facing vertical natural-type surfaces
184	REAL(wp), DIMENSION(1:surf_lsm_v(1)%ns) :: s_flux_lsm_v_south !< flux at south-facing vertical natural-type surfaces
185	REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) :: s_flux_lsm_v_east !< flux at east-facing vertical natural-type surfaces
186	REAL(wp), DIMENSION(1:surf_lsm_v(3)%ns) :: s_flux_lsm_v_west !< flux at west-facing vertical natural-type surfaces
187	REAL(wp), DIMENSION(1:surf_usm_h%ns) :: s_flux_usm_h_up !< flux at horizontal upward-facing urban-type surfaces
188	REAL(wp), DIMENSION(1:surf_usm_v(0)%ns) :: s_flux_usm_v_north !< flux at north-facing vertical urban-type surfaces
189	REAL(wp), DIMENSION(1:surf_usm_v(1)%ns) :: s_flux_usm_v_south !< flux at south-facing vertical urban-type surfaces
190	REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) :: s_flux_usm_v_east !< flux at east-facing vertical urban-type surfaces
191	REAL(wp), DIMENSION(1:surf_usm_v(3)%ns) :: s_flux_usm_v_west !< flux at west-facing vertical urban-type surfaces
192	REAL(wp), DIMENSION(1:surf_def_h(2)%ns) :: s_flux_t !< flux at model top
193	#if defined( __nopointer )
194	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: s !< treated scalar
195	#else
196	REAL(wp), DIMENSION(:,:,:), POINTER :: s !< treated scalar
197	#endif
198
199	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, m) &
200	!$ACC PRIVATE(surf_e, surf_s, flag, mask_top, mask_bottom) &
201	!$ACC PRIVATE(mask_north, mask_south, mask_west, mask_east) &
202	!$ACC PRESENT(wall_flags_0, kh) &
203	!$ACC PRESENT(s) &
204	!$ACC PRESENT(ddzu, ddzw, drho_air, rho_air_zw) &
205	!$ACC PRESENT(surf_def_h(0:2), surf_def_v(0:3)) &
206	!$ACC PRESENT(surf_lsm_h, surf_lsm_v(0:3)) &
207	!$ACC PRESENT(surf_usm_h, surf_usm_v(0:3)) &
208	!$ACC PRESENT(s_flux_def_h_up, s_flux_def_h_down) &
209	!$ACC PRESENT(s_flux_t) &
210	!$ACC PRESENT(s_flux_def_v_north, s_flux_def_v_south) &
211	!$ACC PRESENT(s_flux_def_v_east, s_flux_def_v_west) &
212	!$ACC PRESENT(s_flux_lsm_h_up) &
213	!$ACC PRESENT(s_flux_lsm_v_north, s_flux_lsm_v_south) &
214	!$ACC PRESENT(s_flux_lsm_v_east, s_flux_lsm_v_west) &
215	!$ACC PRESENT(s_flux_usm_h_up) &
216	!$ACC PRESENT(s_flux_usm_v_north, s_flux_usm_v_south) &
217	!$ACC PRESENT(s_flux_usm_v_east, s_flux_usm_v_west) &
218	!$ACC PRESENT(tend)
219	DO i = nxl, nxr
220	DO j = nys,nyn
221	!
222	!-- Compute horizontal diffusion
223	DO k = nzb+1, nzt
224	!
225	!-- Predetermine flag to mask topography and wall-bounded grid points
226	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
227	!
228	!-- Predetermine flag to mask wall-bounded grid points, equivalent to
229	!-- former s_outer array
230	mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i-1), 0 ) )
231	mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i+1), 0 ) )
232	mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j-1,i), 0 ) )
233	mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j+1,i), 0 ) )
234
235	tend(k,j,i) = tend(k,j,i) &
236	+ 0.5_wp * ( &
237	mask_east * ( kh(k,j,i) + kh(k,j,i+1) ) &
238	* ( s(k,j,i+1) - s(k,j,i) ) &
239	- mask_west * ( kh(k,j,i) + kh(k,j,i-1) ) &
240	* ( s(k,j,i) - s(k,j,i-1) ) &
241	) * ddx2 * flag &
242	+ 0.5_wp * ( &
243	mask_north * ( kh(k,j,i) + kh(k,j+1,i) ) &
244	* ( s(k,j+1,i) - s(k,j,i) ) &
245	- mask_south * ( kh(k,j,i) + kh(k,j-1,i) ) &
246	* ( s(k,j,i) - s(k,j-1,i) ) &
247	) * ddy2 * flag
248	ENDDO
249
250	!
251	!-- Apply prescribed horizontal wall heatflux where necessary. First,
252	!-- determine start and end index for respective (j,i)-index. Please
253	!-- note, in the flat case following loop will not be entered, as
254	!-- surf_s=1 and surf_e=0. Furtermore, note, no vertical natural surfaces
255	!-- so far.
256	!-- First, for default-type surfaces
257	!-- North-facing vertical default-type surfaces
258	surf_s = surf_def_v(0)%start_index(j,i)
259	surf_e = surf_def_v(0)%end_index(j,i)
260	DO m = surf_s, surf_e
261	k = surf_def_v(0)%k(m)
262	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_north(m) * ddy
263	ENDDO
264	!
265	!-- South-facing vertical default-type surfaces
266	surf_s = surf_def_v(1)%start_index(j,i)
267	surf_e = surf_def_v(1)%end_index(j,i)
268	DO m = surf_s, surf_e
269	k = surf_def_v(1)%k(m)
270	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_south(m) * ddy
271	ENDDO
272	!
273	!-- East-facing vertical default-type surfaces
274	surf_s = surf_def_v(2)%start_index(j,i)
275	surf_e = surf_def_v(2)%end_index(j,i)
276	DO m = surf_s, surf_e
277	k = surf_def_v(2)%k(m)
278	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_east(m) * ddx
279	ENDDO
280	!
281	!-- West-facing vertical default-type surfaces
282	surf_s = surf_def_v(3)%start_index(j,i)
283	surf_e = surf_def_v(3)%end_index(j,i)
284	DO m = surf_s, surf_e
285	k = surf_def_v(3)%k(m)
286	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_west(m) * ddx
287	ENDDO
288	!
289	!-- Now, for natural-type surfaces.
290	!-- North-facing
291	surf_s = surf_lsm_v(0)%start_index(j,i)
292	surf_e = surf_lsm_v(0)%end_index(j,i)
293	DO m = surf_s, surf_e
294	k = surf_lsm_v(0)%k(m)
295	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_north(m) * ddy
296	ENDDO
297	!
298	!-- South-facing
299	surf_s = surf_lsm_v(1)%start_index(j,i)
300	surf_e = surf_lsm_v(1)%end_index(j,i)
301	DO m = surf_s, surf_e
302	k = surf_lsm_v(1)%k(m)
303	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_south(m) * ddy
304	ENDDO
305	!
306	!-- East-facing
307	surf_s = surf_lsm_v(2)%start_index(j,i)
308	surf_e = surf_lsm_v(2)%end_index(j,i)
309	DO m = surf_s, surf_e
310	k = surf_lsm_v(2)%k(m)
311	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_east(m) * ddx
312	ENDDO
313	!
314	!-- West-facing
315	surf_s = surf_lsm_v(3)%start_index(j,i)
316	surf_e = surf_lsm_v(3)%end_index(j,i)
317	DO m = surf_s, surf_e
318	k = surf_lsm_v(3)%k(m)
319	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_west(m) * ddx
320	ENDDO
321	!
322	!-- Now, for urban-type surfaces.
323	!-- North-facing
324	surf_s = surf_usm_v(0)%start_index(j,i)
325	surf_e = surf_usm_v(0)%end_index(j,i)
326	DO m = surf_s, surf_e
327	k = surf_usm_v(0)%k(m)
328	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_north(m) * ddy
329	ENDDO
330	!
331	!-- South-facing
332	surf_s = surf_usm_v(1)%start_index(j,i)
333	surf_e = surf_usm_v(1)%end_index(j,i)
334	DO m = surf_s, surf_e
335	k = surf_usm_v(1)%k(m)
336	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_south(m) * ddy
337	ENDDO
338	!
339	!-- East-facing
340	surf_s = surf_usm_v(2)%start_index(j,i)
341	surf_e = surf_usm_v(2)%end_index(j,i)
342	DO m = surf_s, surf_e
343	k = surf_usm_v(2)%k(m)
344	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_east(m) * ddx
345	ENDDO
346	!
347	!-- West-facing
348	surf_s = surf_usm_v(3)%start_index(j,i)
349	surf_e = surf_usm_v(3)%end_index(j,i)
350	DO m = surf_s, surf_e
351	k = surf_usm_v(3)%k(m)
352	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_west(m) * ddx
353	ENDDO
354
355	!
356	!-- Compute vertical diffusion. In case that surface fluxes have been
357	!-- prescribed or computed at bottom and/or top, index k starts/ends at
358	!-- nzb+2 or nzt-1, respectively. Model top is also mask if top flux
359	!-- is given.
360	DO k = nzb+1, nzt
361	!
362	!-- Determine flags to mask topography below and above. Flag 0 is
363	!-- used to mask topography in general, and flag 8 implies
364	!-- information about use_surface_fluxes. Flag 9 is used to control
365	!-- flux at model top.
366	mask_bottom = MERGE( 1.0_wp, 0.0_wp, &
367	BTEST( wall_flags_0(k-1,j,i), 8 ) )
368	mask_top = MERGE( 1.0_wp, 0.0_wp, &
369	BTEST( wall_flags_0(k+1,j,i), 8 ) ) * &
370	MERGE( 1.0_wp, 0.0_wp, &
371	BTEST( wall_flags_0(k+1,j,i), 9 ) )
372	flag = MERGE( 1.0_wp, 0.0_wp, &
373	BTEST( wall_flags_0(k,j,i), 0 ) )
374
375	tend(k,j,i) = tend(k,j,i) &
376	+ 0.5_wp * ( &
377	( kh(k,j,i) + kh(k+1,j,i) ) * &
378	( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
379	* rho_air_zw(k) &
380	* mask_top &
381	- ( kh(k,j,i) + kh(k-1,j,i) ) * &
382	( s(k,j,i)-s(k-1,j,i) ) * ddzu(k) &
383	* rho_air_zw(k-1) &
384	* mask_bottom &
385	) * ddzw(k) * drho_air(k) &
386	* flag
387	ENDDO
388
389	!
390	!-- Vertical diffusion at horizontal walls.
391	IF ( use_surface_fluxes ) THEN
392	!
393	!-- Default-type surfaces, upward-facing
394	surf_s = surf_def_h(0)%start_index(j,i)
395	surf_e = surf_def_h(0)%end_index(j,i)
396	DO m = surf_s, surf_e
397
398	k = surf_def_h(0)%k(m)
399	tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m) &
400	* ddzw(k) * drho_air(k)
401
402	ENDDO
403	!
404	!-- Default-type surfaces, downward-facing
405	surf_s = surf_def_h(1)%start_index(j,i)
406	surf_e = surf_def_h(1)%end_index(j,i)
407	DO m = surf_s, surf_e
408
409	k = surf_def_h(1)%k(m)
410	tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m) &
411	* ddzw(k) * drho_air(k)
412
413	ENDDO
414	!
415	!-- Natural-type surfaces, upward-facing
416	surf_s = surf_lsm_h%start_index(j,i)
417	surf_e = surf_lsm_h%end_index(j,i)
418	DO m = surf_s, surf_e
419
420	k = surf_lsm_h%k(m)
421	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_up(m) &
422	* ddzw(k) * drho_air(k)
423
424	ENDDO
425	!
426	!-- Urban-type surfaces, upward-facing
427	surf_s = surf_usm_h%start_index(j,i)
428	surf_e = surf_usm_h%end_index(j,i)
429	DO m = surf_s, surf_e
430
431	k = surf_usm_h%k(m)
432	tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_up(m) &
433	* ddzw(k) * drho_air(k)
434
435	ENDDO
436
437	ENDIF
438	!
439	!-- Vertical diffusion at the last computational gridpoint along z-direction
440	IF ( use_top_fluxes ) THEN
441	surf_s = surf_def_h(2)%start_index(j,i)
442	surf_e = surf_def_h(2)%end_index(j,i)
443	DO m = surf_s, surf_e
444
445	k = surf_def_h(2)%k(m)
446	tend(k,j,i) = tend(k,j,i) &
447	+ ( - s_flux_t(m) ) * ddzw(k) * drho_air(k)
448	ENDDO
449	ENDIF
450
451	ENDDO
452	ENDDO
453
454	END SUBROUTINE diffusion_s
455
456	!------------------------------------------------------------------------------!
457	! Description:
458	! ------------
459	!> Call for grid point i,j
460	!------------------------------------------------------------------------------!
461	SUBROUTINE diffusion_s_ij( i, j, s, &
462	s_flux_def_h_up, s_flux_def_h_down, &
463	s_flux_t, &
464	s_flux_lsm_h_up, s_flux_usm_h_up, &
465	s_flux_def_v_north, s_flux_def_v_south, &
466	s_flux_def_v_east, s_flux_def_v_west, &
467	s_flux_lsm_v_north, s_flux_lsm_v_south, &
468	s_flux_lsm_v_east, s_flux_lsm_v_west, &
469	s_flux_usm_v_north, s_flux_usm_v_south, &
470	s_flux_usm_v_east, s_flux_usm_v_west )
471
472	USE arrays_3d, &
473	ONLY: ddzu, ddzw, kh, tend, drho_air, rho_air_zw
474
475	USE control_parameters, &
476	ONLY: use_surface_fluxes, use_top_fluxes
477
478	USE grid_variables, &
479	ONLY: ddx, ddx2, ddy, ddy2
480
481	USE indices, &
482	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
483
484	USE kinds
485
486	USE surface_mod, &
487	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
488	surf_usm_v
489
490	IMPLICIT NONE
491
492	INTEGER(iwp) :: i !< running index x direction
493	INTEGER(iwp) :: j !< running index y direction
494	INTEGER(iwp) :: k !< running index z direction
495	INTEGER(iwp) :: m !< running index surface elements
496	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
497	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
498
499	REAL(wp) :: flag !< flag to mask topography grid points
500	REAL(wp) :: mask_bottom !< flag to mask vertical upward-facing surface
501	REAL(wp) :: mask_east !< flag to mask vertical surface east of the grid point
502	REAL(wp) :: mask_north !< flag to mask vertical surface north of the grid point
503	REAL(wp) :: mask_south !< flag to mask vertical surface south of the grid point
504	REAL(wp) :: mask_west !< flag to mask vertical surface west of the grid point
505	REAL(wp) :: mask_top !< flag to mask vertical downward-facing surface
506
507	REAL(wp), DIMENSION(1:surf_def_v(0)%ns) :: s_flux_def_v_north !< flux at north-facing vertical default-type surfaces
508	REAL(wp), DIMENSION(1:surf_def_v(1)%ns) :: s_flux_def_v_south !< flux at south-facing vertical default-type surfaces
509	REAL(wp), DIMENSION(1:surf_def_v(2)%ns) :: s_flux_def_v_east !< flux at east-facing vertical default-type surfaces
510	REAL(wp), DIMENSION(1:surf_def_v(3)%ns) :: s_flux_def_v_west !< flux at west-facing vertical default-type surfaces
511	REAL(wp), DIMENSION(1:surf_def_h(0)%ns) :: s_flux_def_h_up !< flux at horizontal upward-facing default-type surfaces
512	REAL(wp), DIMENSION(1:surf_def_h(1)%ns) :: s_flux_def_h_down !< flux at horizontal donwward-facing default-type surfaces
513	REAL(wp), DIMENSION(1:surf_lsm_h%ns) :: s_flux_lsm_h_up !< flux at horizontal upward-facing natural-type surfaces
514	REAL(wp), DIMENSION(1:surf_lsm_v(0)%ns) :: s_flux_lsm_v_north !< flux at north-facing vertical urban-type surfaces
515	REAL(wp), DIMENSION(1:surf_lsm_v(1)%ns) :: s_flux_lsm_v_south !< flux at south-facing vertical urban-type surfaces
516	REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) :: s_flux_lsm_v_east !< flux at east-facing vertical urban-type surfaces
517	REAL(wp), DIMENSION(1:surf_lsm_v(3)%ns) :: s_flux_lsm_v_west !< flux at west-facing vertical urban-type surfaces
518	REAL(wp), DIMENSION(1:surf_usm_h%ns) :: s_flux_usm_h_up !< flux at horizontal upward-facing urban-type surfaces
519	REAL(wp), DIMENSION(1:surf_usm_v(0)%ns) :: s_flux_usm_v_north !< flux at north-facing vertical urban-type surfaces
520	REAL(wp), DIMENSION(1:surf_usm_v(1)%ns) :: s_flux_usm_v_south !< flux at south-facing vertical urban-type surfaces
521	REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) :: s_flux_usm_v_east !< flux at east-facing vertical urban-type surfaces
522	REAL(wp), DIMENSION(1:surf_usm_v(3)%ns) :: s_flux_usm_v_west !< flux at west-facing vertical urban-type surfaces
523	REAL(wp), DIMENSION(1:surf_def_h(2)%ns) :: s_flux_t !< flux at model top
524	#if defined( __nopointer )
525	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: s !< treated scalar
526	#else
527	REAL(wp), DIMENSION(:,:,:), POINTER :: s !< treated scalar
528	#endif
529
530	!
531	!-- Compute horizontal diffusion
532	DO k = nzb+1, nzt
533	!
534	!-- Predetermine flag to mask topography and wall-bounded grid points
535	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
536	!
537	!-- Predetermine flag to mask wall-bounded grid points, equivalent to
538	!-- former s_outer array
539	mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i-1), 0 ) )
540	mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i+1), 0 ) )
541	mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j-1,i), 0 ) )
542	mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j+1,i), 0 ) )
543	!
544	!-- Finally, determine flag to mask both topography itself as well
545	!-- as wall-bounded grid points, which will be treated further below
546
547	tend(k,j,i) = tend(k,j,i) &
548	+ 0.5_wp * ( &
549	mask_east * ( kh(k,j,i) + kh(k,j,i+1) ) &
550	* ( s(k,j,i+1) - s(k,j,i) ) &
551	- mask_west * ( kh(k,j,i) + kh(k,j,i-1) ) &
552	* ( s(k,j,i) - s(k,j,i-1) ) &
553	) * ddx2 * flag &
554	+ 0.5_wp * ( &
555	mask_north * ( kh(k,j,i) + kh(k,j+1,i) ) &
556	* ( s(k,j+1,i) - s(k,j,i) ) &
557	- mask_south * ( kh(k,j,i) + kh(k,j-1,i) ) &
558	* ( s(k,j,i) - s(k,j-1,i) ) &
559	) * ddy2 * flag
560	ENDDO
561
562	!
563	!-- Apply prescribed horizontal wall heatflux where necessary. First,
564	!-- determine start and end index for respective (j,i)-index. Please
565	!-- note, in the flat case following loops will not be entered, as
566	!-- surf_s=1 and surf_e=0. Furtermore, note, no vertical natural surfaces
567	!-- so far.
568	!-- First, for default-type surfaces
569	!-- North-facing vertical default-type surfaces
570	surf_s = surf_def_v(0)%start_index(j,i)
571	surf_e = surf_def_v(0)%end_index(j,i)
572	DO m = surf_s, surf_e
573	k = surf_def_v(0)%k(m)
574	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_north(m) * ddy
575	ENDDO
576	!
577	!-- South-facing vertical default-type surfaces
578	surf_s = surf_def_v(1)%start_index(j,i)
579	surf_e = surf_def_v(1)%end_index(j,i)
580	DO m = surf_s, surf_e
581	k = surf_def_v(1)%k(m)
582	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_south(m) * ddy
583	ENDDO
584	!
585	!-- East-facing vertical default-type surfaces
586	surf_s = surf_def_v(2)%start_index(j,i)
587	surf_e = surf_def_v(2)%end_index(j,i)
588	DO m = surf_s, surf_e
589	k = surf_def_v(2)%k(m)
590	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_east(m) * ddx
591	ENDDO
592	!
593	!-- West-facing vertical default-type surfaces
594	surf_s = surf_def_v(3)%start_index(j,i)
595	surf_e = surf_def_v(3)%end_index(j,i)
596	DO m = surf_s, surf_e
597	k = surf_def_v(3)%k(m)
598	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_west(m) * ddx
599	ENDDO
600	!
601	!-- Now, for natural-type surfaces
602	!-- North-facing
603	surf_s = surf_lsm_v(0)%start_index(j,i)
604	surf_e = surf_lsm_v(0)%end_index(j,i)
605	DO m = surf_s, surf_e
606	k = surf_lsm_v(0)%k(m)
607	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_north(m) * ddy
608	ENDDO
609	!
610	!-- South-facing
611	surf_s = surf_lsm_v(1)%start_index(j,i)
612	surf_e = surf_lsm_v(1)%end_index(j,i)
613	DO m = surf_s, surf_e
614	k = surf_lsm_v(1)%k(m)
615	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_south(m) * ddy
616	ENDDO
617	!
618	!-- East-facing
619	surf_s = surf_lsm_v(2)%start_index(j,i)
620	surf_e = surf_lsm_v(2)%end_index(j,i)
621	DO m = surf_s, surf_e
622	k = surf_lsm_v(2)%k(m)
623	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_east(m) * ddx
624	ENDDO
625	!
626	!-- West-facing
627	surf_s = surf_lsm_v(3)%start_index(j,i)
628	surf_e = surf_lsm_v(3)%end_index(j,i)
629	DO m = surf_s, surf_e
630	k = surf_lsm_v(3)%k(m)
631	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_west(m) * ddx
632	ENDDO
633	!
634	!-- Now, for urban-type surfaces
635	!-- North-facing
636	surf_s = surf_usm_v(0)%start_index(j,i)
637	surf_e = surf_usm_v(0)%end_index(j,i)
638	DO m = surf_s, surf_e
639	k = surf_usm_v(0)%k(m)
640	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_north(m) * ddy
641	ENDDO
642	!
643	!-- South-facing
644	surf_s = surf_usm_v(1)%start_index(j,i)
645	surf_e = surf_usm_v(1)%end_index(j,i)
646	DO m = surf_s, surf_e
647	k = surf_usm_v(1)%k(m)
648	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_south(m) * ddy
649	ENDDO
650	!
651	!-- East-facing
652	surf_s = surf_usm_v(2)%start_index(j,i)
653	surf_e = surf_usm_v(2)%end_index(j,i)
654	DO m = surf_s, surf_e
655	k = surf_usm_v(2)%k(m)
656	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_east(m) * ddx
657	ENDDO
658	!
659	!-- West-facing
660	surf_s = surf_usm_v(3)%start_index(j,i)
661	surf_e = surf_usm_v(3)%end_index(j,i)
662	DO m = surf_s, surf_e
663	k = surf_usm_v(3)%k(m)
664	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_west(m) * ddx
665	ENDDO
666
667
668	!
669	!-- Compute vertical diffusion. In case that surface fluxes have been
670	!-- prescribed or computed at bottom and/or top, index k starts/ends at
671	!-- nzb+2 or nzt-1, respectively. Model top is also mask if top flux
672	!-- is given.
673	DO k = nzb+1, nzt
674	!
675	!-- Determine flags to mask topography below and above. Flag 0 is
676	!-- used to mask topography in general, and flag 8 implies
677	!-- information about use_surface_fluxes. Flag 9 is used to control
678	!-- flux at model top.
679	mask_bottom = MERGE( 1.0_wp, 0.0_wp, &
680	BTEST( wall_flags_0(k-1,j,i), 8 ) )
681	mask_top = MERGE( 1.0_wp, 0.0_wp, &
682	BTEST( wall_flags_0(k+1,j,i), 8 ) ) * &
683	MERGE( 1.0_wp, 0.0_wp, &
684	BTEST( wall_flags_0(k+1,j,i), 9 ) )
685	flag = MERGE( 1.0_wp, 0.0_wp, &
686	BTEST( wall_flags_0(k,j,i), 0 ) )
687
688	tend(k,j,i) = tend(k,j,i) &
689	+ 0.5_wp * ( &
690	( kh(k,j,i) + kh(k+1,j,i) ) * &
691	( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
692	* rho_air_zw(k) &
693	* mask_top &
694	- ( kh(k,j,i) + kh(k-1,j,i) ) * &
695	( s(k,j,i)-s(k-1,j,i) ) * ddzu(k) &
696	* rho_air_zw(k-1) &
697	* mask_bottom &
698	) * ddzw(k) * drho_air(k) &
699	* flag
700	ENDDO
701
702	!
703	!-- Vertical diffusion at horizontal walls.
704	!-- TO DO: Adjust for downward facing walls and mask already in main loop
705	IF ( use_surface_fluxes ) THEN
706	!
707	!-- Default-type surfaces, upward-facing
708	surf_s = surf_def_h(0)%start_index(j,i)
709	surf_e = surf_def_h(0)%end_index(j,i)
710	DO m = surf_s, surf_e
711
712	k = surf_def_h(0)%k(m)
713
714	tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m) &
715	* ddzw(k) * drho_air(k)
716	ENDDO
717	!
718	!-- Default-type surfaces, downward-facing
719	surf_s = surf_def_h(1)%start_index(j,i)
720	surf_e = surf_def_h(1)%end_index(j,i)
721	DO m = surf_s, surf_e
722
723	k = surf_def_h(1)%k(m)
724
725	tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m) &
726	* ddzw(k) * drho_air(k)
727	ENDDO
728	!
729	!-- Natural-type surfaces, upward-facing
730	surf_s = surf_lsm_h%start_index(j,i)
731	surf_e = surf_lsm_h%end_index(j,i)
732	DO m = surf_s, surf_e
733	k = surf_lsm_h%k(m)
734
735	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_up(m) &
736	* ddzw(k) * drho_air(k)
737	ENDDO
738	!
739	!-- Urban-type surfaces, upward-facing
740	surf_s = surf_usm_h%start_index(j,i)
741	surf_e = surf_usm_h%end_index(j,i)
742	DO m = surf_s, surf_e
743	k = surf_usm_h%k(m)
744
745	tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_up(m) &
746	* ddzw(k) * drho_air(k)
747	ENDDO
748	ENDIF
749	!
750	!-- Vertical diffusion at the last computational gridpoint along z-direction
751	IF ( use_top_fluxes ) THEN
752	surf_s = surf_def_h(2)%start_index(j,i)
753	surf_e = surf_def_h(2)%end_index(j,i)
754	DO m = surf_s, surf_e
755
756	k = surf_def_h(2)%k(m)
757	tend(k,j,i) = tend(k,j,i) &
758	+ ( - s_flux_t(m) ) * ddzw(k) * drho_air(k)
759	ENDDO
760	ENDIF
761
762	END SUBROUTINE diffusion_s_ij
763
764	END MODULE diffusion_s_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |