Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |