Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |