Home

Context Navigation

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

Last change on this file since 2232 was 2232, checked in by suehring, 7 years ago
Adjustments according new topography and surface-modelling concept implemented
Property svn:keywords set to `Id`
File size: 33.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |