Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |