Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |