1 | MODULE advec_ws |
---|
2 | |
---|
3 | !-----------------------------------------------------------------------------! |
---|
4 | ! Current revisions: |
---|
5 | ! ------------------ |
---|
6 | ! ws-scheme also work with topography in combination with vector version. |
---|
7 | ! ws-scheme also work with outflow boundaries in combination with |
---|
8 | ! vector version. |
---|
9 | |
---|
10 | ! Degradation of the applied order of scheme is now steered by multiplying with |
---|
11 | ! Integer wall_flags_0. 2nd order scheme, WS3 and WS5 are calculated on each |
---|
12 | ! grid point and mulitplied with the appropriate flag. |
---|
13 | ! 2nd order numerical dissipation term changed. Now the appropriate 2nd order |
---|
14 | ! term derived according to the 4th and 6th order terms is applied. It turns |
---|
15 | ! out that diss_2nd does not provide sufficient dissipation near walls. |
---|
16 | ! Therefore, the function diss_2nd is removed. |
---|
17 | |
---|
18 | ! Near walls a divergence correction is necessary to overcome numerical |
---|
19 | ! instabilities due to too less divergence reduction of the velocity field. |
---|
20 | ! |
---|
21 | ! boundary_flags and logicals steering the degradation are removed. |
---|
22 | ! Empty SUBROUTINE local_diss removed. |
---|
23 | ! |
---|
24 | ! Further formatting adjustments. |
---|
25 | ! |
---|
26 | ! Former revisions: |
---|
27 | ! ----------------- |
---|
28 | ! $Id: advec_ws.f90 862 2012-03-26 14:21:38Z suehring $ |
---|
29 | ! |
---|
30 | ! 801 2012-01-10 17:30:36Z suehring |
---|
31 | ! Bugfix concerning OpenMP parallelization. Summation of sums_wsus_ws_l, |
---|
32 | ! sums_wsvs_ws_l, sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wspts_ws_l, |
---|
33 | ! sums_wsqs_ws_l, sums_wssas_ws_l is now thread-safe by adding an additional |
---|
34 | ! dimension. |
---|
35 | ! |
---|
36 | ! 743 2011-08-18 16:10:16Z suehring |
---|
37 | ! Evaluation of turbulent fluxes with WS-scheme only for the whole model |
---|
38 | ! domain. Therefor dimension of arrays needed for statistical evaluation |
---|
39 | ! decreased. |
---|
40 | ! |
---|
41 | ! 736 2011-08-17 14:13:26Z suehring |
---|
42 | ! Bugfix concerning OpenMP parallelization. i_omp introduced, because first |
---|
43 | ! index where fluxes on left side have to be calculated explicitly is |
---|
44 | ! different on each thread. Furthermore the swapping of fluxes is now |
---|
45 | ! thread-safe by adding an additional dimension. |
---|
46 | ! |
---|
47 | ! 713 2011-03-30 14:21:21Z suehring |
---|
48 | ! File reformatted. |
---|
49 | ! Bugfix in vertical advection of w concerning the optimized version for |
---|
50 | ! vector architecture. |
---|
51 | ! Constants adv_mom_3, adv_mom_5, adv_sca_5, adv_sca_3 reformulated as |
---|
52 | ! broken numbers. |
---|
53 | ! |
---|
54 | ! 709 2011-03-30 09:31:40Z raasch |
---|
55 | ! formatting adjustments |
---|
56 | ! |
---|
57 | ! 705 2011-03-25 11:21:43 Z suehring $ |
---|
58 | ! Bugfix in declaration of logicals concerning outflow boundaries. |
---|
59 | ! |
---|
60 | ! 411 2009-12-11 12:31:43 Z suehring |
---|
61 | ! Allocation of weight_substep moved to init_3d_model. |
---|
62 | ! Declaration of ws_scheme_sca and ws_scheme_mom moved to check_parameters. |
---|
63 | ! Setting bc for the horizontal velocity variances added (moved from |
---|
64 | ! flow_statistics). |
---|
65 | ! |
---|
66 | ! Initial revision |
---|
67 | ! |
---|
68 | ! 411 2009-12-11 12:31:43 Z suehring |
---|
69 | ! |
---|
70 | ! Description: |
---|
71 | ! ------------ |
---|
72 | ! Advection scheme for scalars and momentum using the flux formulation of |
---|
73 | ! Wicker and Skamarock 5th order. Additionally the module contains of a |
---|
74 | ! routine using for initialisation and steering of the statical evaluation. |
---|
75 | ! The computation of turbulent fluxes takes place inside the advection |
---|
76 | ! routines. |
---|
77 | ! Near non-cyclic boundaries the order of the applied advection scheme is |
---|
78 | ! degraded. |
---|
79 | ! A divergence correction is applied. It is necessary for topography, since |
---|
80 | ! the divergence is not sufficiently reduced, resulting in erroneous fluxes and |
---|
81 | ! partly numerical instabilities. |
---|
82 | !-----------------------------------------------------------------------------! |
---|
83 | |
---|
84 | PRIVATE |
---|
85 | PUBLIC advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws, & |
---|
86 | ws_init, ws_statistics |
---|
87 | |
---|
88 | INTERFACE ws_init |
---|
89 | MODULE PROCEDURE ws_init |
---|
90 | END INTERFACE ws_init |
---|
91 | |
---|
92 | INTERFACE ws_statistics |
---|
93 | MODULE PROCEDURE ws_statistics |
---|
94 | END INTERFACE ws_statistics |
---|
95 | |
---|
96 | INTERFACE advec_s_ws |
---|
97 | MODULE PROCEDURE advec_s_ws |
---|
98 | MODULE PROCEDURE advec_s_ws_ij |
---|
99 | END INTERFACE advec_s_ws |
---|
100 | |
---|
101 | INTERFACE advec_u_ws |
---|
102 | MODULE PROCEDURE advec_u_ws |
---|
103 | MODULE PROCEDURE advec_u_ws_ij |
---|
104 | END INTERFACE advec_u_ws |
---|
105 | |
---|
106 | INTERFACE advec_v_ws |
---|
107 | MODULE PROCEDURE advec_v_ws |
---|
108 | MODULE PROCEDURE advec_v_ws_ij |
---|
109 | END INTERFACE advec_v_ws |
---|
110 | |
---|
111 | INTERFACE advec_w_ws |
---|
112 | MODULE PROCEDURE advec_w_ws |
---|
113 | MODULE PROCEDURE advec_w_ws_ij |
---|
114 | END INTERFACE advec_w_ws |
---|
115 | |
---|
116 | CONTAINS |
---|
117 | |
---|
118 | |
---|
119 | !------------------------------------------------------------------------------! |
---|
120 | ! Initialization of WS-scheme |
---|
121 | !------------------------------------------------------------------------------! |
---|
122 | SUBROUTINE ws_init |
---|
123 | |
---|
124 | USE arrays_3d |
---|
125 | USE constants |
---|
126 | USE control_parameters |
---|
127 | USE indices |
---|
128 | USE pegrid |
---|
129 | USE statistics |
---|
130 | |
---|
131 | ! |
---|
132 | !-- Set the appropriate factors for scalar and momentum advection. |
---|
133 | adv_sca_5 = 1./60. |
---|
134 | adv_sca_3 = 1./12. |
---|
135 | adv_sca_1 = 1./2. |
---|
136 | adv_mom_5 = 1./120. |
---|
137 | adv_mom_3 = 1./24. |
---|
138 | adv_mom_1 = 1./4. |
---|
139 | ! |
---|
140 | !-- Arrays needed for statical evaluation of fluxes. |
---|
141 | IF ( ws_scheme_mom ) THEN |
---|
142 | |
---|
143 | ALLOCATE( sums_wsus_ws_l(nzb:nzt+1,0:threads_per_task-1), & |
---|
144 | sums_wsvs_ws_l(nzb:nzt+1,0:threads_per_task-1), & |
---|
145 | sums_us2_ws_l(nzb:nzt+1,0:threads_per_task-1), & |
---|
146 | sums_vs2_ws_l(nzb:nzt+1,0:threads_per_task-1), & |
---|
147 | sums_ws2_ws_l(nzb:nzt+1,0:threads_per_task-1) ) |
---|
148 | |
---|
149 | sums_wsus_ws_l = 0.0 |
---|
150 | sums_wsvs_ws_l = 0.0 |
---|
151 | sums_us2_ws_l = 0.0 |
---|
152 | sums_vs2_ws_l = 0.0 |
---|
153 | sums_ws2_ws_l = 0.0 |
---|
154 | |
---|
155 | ENDIF |
---|
156 | |
---|
157 | IF ( ws_scheme_sca ) THEN |
---|
158 | |
---|
159 | ALLOCATE( sums_wspts_ws_l(nzb:nzt+1,0:threads_per_task-1) ) |
---|
160 | sums_wspts_ws_l = 0.0 |
---|
161 | |
---|
162 | IF ( humidity .OR. passive_scalar ) THEN |
---|
163 | ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1,0:threads_per_task-1) ) |
---|
164 | sums_wsqs_ws_l = 0.0 |
---|
165 | ENDIF |
---|
166 | |
---|
167 | IF ( ocean ) THEN |
---|
168 | ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:threads_per_task-1) ) |
---|
169 | sums_wssas_ws_l = 0.0 |
---|
170 | ENDIF |
---|
171 | |
---|
172 | ENDIF |
---|
173 | |
---|
174 | ! |
---|
175 | !-- Arrays needed for reasons of speed optimization for cache and noopt |
---|
176 | !-- version. For the vector version the buffer arrays are not necessary, |
---|
177 | !-- because the the fluxes can swapped directly inside the loops of the |
---|
178 | !-- advection routines. |
---|
179 | IF ( loop_optimization /= 'vector' ) THEN |
---|
180 | |
---|
181 | IF ( ws_scheme_mom ) THEN |
---|
182 | |
---|
183 | ALLOCATE( flux_s_u(nzb+1:nzt,0:threads_per_task-1), & |
---|
184 | flux_s_v(nzb+1:nzt,0:threads_per_task-1), & |
---|
185 | flux_s_w(nzb+1:nzt,0:threads_per_task-1), & |
---|
186 | diss_s_u(nzb+1:nzt,0:threads_per_task-1), & |
---|
187 | diss_s_v(nzb+1:nzt,0:threads_per_task-1), & |
---|
188 | diss_s_w(nzb+1:nzt,0:threads_per_task-1) ) |
---|
189 | ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & |
---|
190 | flux_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & |
---|
191 | flux_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & |
---|
192 | diss_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & |
---|
193 | diss_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & |
---|
194 | diss_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) |
---|
195 | |
---|
196 | ENDIF |
---|
197 | |
---|
198 | IF ( ws_scheme_sca ) THEN |
---|
199 | |
---|
200 | ALLOCATE( flux_s_pt(nzb+1:nzt,0:threads_per_task-1), & |
---|
201 | flux_s_e(nzb+1:nzt,0:threads_per_task-1), & |
---|
202 | diss_s_pt(nzb+1:nzt,0:threads_per_task-1), & |
---|
203 | diss_s_e(nzb+1:nzt,0:threads_per_task-1) ) |
---|
204 | ALLOCATE( flux_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & |
---|
205 | flux_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & |
---|
206 | diss_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & |
---|
207 | diss_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) |
---|
208 | |
---|
209 | IF ( humidity .OR. passive_scalar ) THEN |
---|
210 | ALLOCATE( flux_s_q(nzb+1:nzt,0:threads_per_task-1), & |
---|
211 | diss_s_q(nzb+1:nzt,0:threads_per_task-1) ) |
---|
212 | ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & |
---|
213 | diss_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) |
---|
214 | ENDIF |
---|
215 | |
---|
216 | IF ( ocean ) THEN |
---|
217 | ALLOCATE( flux_s_sa(nzb+1:nzt,0:threads_per_task-1), & |
---|
218 | diss_s_sa(nzb+1:nzt,0:threads_per_task-1) ) |
---|
219 | ALLOCATE( flux_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & |
---|
220 | diss_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) |
---|
221 | ENDIF |
---|
222 | |
---|
223 | ENDIF |
---|
224 | |
---|
225 | ENDIF |
---|
226 | |
---|
227 | END SUBROUTINE ws_init |
---|
228 | |
---|
229 | |
---|
230 | !------------------------------------------------------------------------------! |
---|
231 | ! Initialize variables used for storing statistic qauntities (fluxes, variances) |
---|
232 | !------------------------------------------------------------------------------! |
---|
233 | SUBROUTINE ws_statistics |
---|
234 | |
---|
235 | USE control_parameters |
---|
236 | USE statistics |
---|
237 | |
---|
238 | IMPLICIT NONE |
---|
239 | |
---|
240 | ! |
---|
241 | !-- The arrays needed for statistical evaluation are set to to 0 at the |
---|
242 | !-- beginning of prognostic_equations. |
---|
243 | IF ( ws_scheme_mom ) THEN |
---|
244 | sums_wsus_ws_l = 0.0 |
---|
245 | sums_wsvs_ws_l = 0.0 |
---|
246 | sums_us2_ws_l = 0.0 |
---|
247 | sums_vs2_ws_l = 0.0 |
---|
248 | sums_ws2_ws_l = 0.0 |
---|
249 | ENDIF |
---|
250 | |
---|
251 | IF ( ws_scheme_sca ) THEN |
---|
252 | sums_wspts_ws_l = 0.0 |
---|
253 | IF ( humidity .OR. passive_scalar ) sums_wsqs_ws_l = 0.0 |
---|
254 | IF ( ocean ) sums_wssas_ws_l = 0.0 |
---|
255 | |
---|
256 | ENDIF |
---|
257 | |
---|
258 | END SUBROUTINE ws_statistics |
---|
259 | |
---|
260 | |
---|
261 | !------------------------------------------------------------------------------! |
---|
262 | ! Scalar advection - Call for grid point i,j |
---|
263 | !------------------------------------------------------------------------------! |
---|
264 | SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char,swap_flux_y_local, & |
---|
265 | swap_diss_y_local, swap_flux_x_local, & |
---|
266 | swap_diss_x_local, i_omp, tn ) |
---|
267 | |
---|
268 | USE arrays_3d |
---|
269 | USE constants |
---|
270 | USE control_parameters |
---|
271 | USE grid_variables |
---|
272 | USE indices |
---|
273 | USE pegrid |
---|
274 | USE statistics |
---|
275 | |
---|
276 | IMPLICIT NONE |
---|
277 | |
---|
278 | INTEGER :: i, i_omp, j, k, k_mm, k_pp, k_ppp, tn |
---|
279 | REAL :: diss_d, div, flux_d, u_comp, v_comp |
---|
280 | REAL, DIMENSION(:,:,:), POINTER :: sk |
---|
281 | REAL, DIMENSION(nzb:nzt+1) :: diss_n, diss_r, diss_t, flux_n, & |
---|
282 | flux_r, flux_t |
---|
283 | REAL, DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_diss_y_local, & |
---|
284 | swap_flux_y_local |
---|
285 | REAL, DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: & |
---|
286 | swap_diss_x_local, & |
---|
287 | swap_flux_x_local |
---|
288 | CHARACTER (LEN = *), INTENT(IN) :: sk_char |
---|
289 | |
---|
290 | ! |
---|
291 | !-- Compute southside fluxes of the respective PE bounds. |
---|
292 | IF ( j == nys ) THEN |
---|
293 | ! |
---|
294 | !-- Up to the top of the highest topography. |
---|
295 | DO k = nzb+1, nzb_max |
---|
296 | |
---|
297 | v_comp = v(k,j,i) - v_gtrans |
---|
298 | swap_flux_y_local(k,tn) = v_comp * ( & |
---|
299 | ( 37.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 & |
---|
300 | + 7.0 * IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3 & |
---|
301 | + IBITS(wall_flags_0(k,j,i),3,1) * adv_sca_1 & |
---|
302 | ) * & |
---|
303 | ( sk(k,j,i) + sk(k,j-1,i) ) & |
---|
304 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 & |
---|
305 | + IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3 & |
---|
306 | ) * & |
---|
307 | ( sk(k,j+1,i) + sk(k,j-2,i) ) & |
---|
308 | + ( IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 & |
---|
309 | ) * & |
---|
310 | ( sk(k,j+2,i) + sk(k,j-3,i) ) & |
---|
311 | ) |
---|
312 | |
---|
313 | swap_diss_y_local(k,tn) = -ABS( v_comp ) * ( & |
---|
314 | ( 10.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 & |
---|
315 | + 3.0 * IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3 & |
---|
316 | + IBITS(wall_flags_0(k,j,i),3,1) * adv_sca_1 & |
---|
317 | ) * & |
---|
318 | ( sk(k,j,i) - sk(k,j-1,i) ) & |
---|
319 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 & |
---|
320 | + IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3 & |
---|
321 | ) * & |
---|
322 | ( sk(k,j+1,i) - sk(k,j-2,i) ) & |
---|
323 | + ( IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 & |
---|
324 | ) * & |
---|
325 | ( sk(k,j+2,i) - sk(k,j-3,i) ) & |
---|
326 | ) |
---|
327 | |
---|
328 | ENDDO |
---|
329 | ! |
---|
330 | !-- Above to the top of the highest topography. No degradation necessary. |
---|
331 | DO k = nzb_max+1, nzt |
---|
332 | |
---|
333 | v_comp = v(k,j,i) - v_gtrans |
---|
334 | swap_flux_y_local(k,tn) = v_comp * ( & |
---|
335 | 37.0 * ( sk(k,j,i) + sk(k,j-1,i) ) & |
---|
336 | - 8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) ) & |
---|
337 | + ( sk(k,j+2,i) + sk(k,j-3,i) ) & |
---|
338 | ) * adv_sca_5 |
---|
339 | swap_diss_y_local(k,tn) = -ABS( v_comp ) * ( & |
---|
340 | 10.0 * ( sk(k,j,i) - sk(k,j-1,i) ) & |
---|
341 | - 5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) ) & |
---|
342 | + sk(k,j+2,i) - sk(k,j-3,i) & |
---|
343 | ) * adv_sca_5 |
---|
344 | |
---|
345 | ENDDO |
---|
346 | |
---|
347 | ENDIF |
---|
348 | ! |
---|
349 | !-- Compute leftside fluxes of the respective PE bounds. |
---|
350 | IF ( i == i_omp ) THEN |
---|
351 | |
---|
352 | DO k = nzb+1, nzb_max |
---|
353 | |
---|
354 | u_comp = u(k,j,i) - u_gtrans |
---|
355 | swap_flux_x_local(k,j,tn) = u_comp * ( & |
---|
356 | ( 37.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 & |
---|
357 | + 7.0 * IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3 & |
---|
358 | + IBITS(wall_flags_0(k,j,i),0,1) * adv_sca_1 & |
---|
359 | ) * & |
---|
360 | ( sk(k,j,i) + sk(k,j,i-1) ) & |
---|
361 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 & |
---|
362 | + IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3 & |
---|
363 | ) * & |
---|
364 | ( sk(k,j,i+1) + sk(k,j,i-2) ) & |
---|
365 | + ( IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 & |
---|
366 | ) * ( sk(k,j,i+2) + sk(k,j,i-3) ) & |
---|
367 | ) |
---|
368 | |
---|
369 | swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * ( & |
---|
370 | ( 10.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 & |
---|
371 | + 3.0 * IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3 & |
---|
372 | + IBITS(wall_flags_0(k,j,i),0,1) * adv_sca_1 & |
---|
373 | ) * & |
---|
374 | ( sk(k,j,i) - sk(k,j,i-1) ) & |
---|
375 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 & |
---|
376 | + IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3 & |
---|
377 | ) * & |
---|
378 | ( sk(k,j,i+1) - sk(k,j,i-2) ) & |
---|
379 | + ( IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 & |
---|
380 | ) * & |
---|
381 | ( sk(k,j,i+2) - sk(k,j,i-3) ) & |
---|
382 | ) |
---|
383 | |
---|
384 | ENDDO |
---|
385 | |
---|
386 | DO k = nzb_max+1, nzt |
---|
387 | |
---|
388 | u_comp = u(k,j,i) - u_gtrans |
---|
389 | swap_flux_x_local(k,j,tn) = u_comp * ( & |
---|
390 | 37.0 * ( sk(k,j,i) + sk(k,j,i-1) ) & |
---|
391 | - 8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) ) & |
---|
392 | + ( sk(k,j,i+2) + sk(k,j,i-3) ) & |
---|
393 | ) * adv_sca_5 |
---|
394 | |
---|
395 | swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * ( & |
---|
396 | 10.0 * ( sk(k,j,i) - sk(k,j,i-1) ) & |
---|
397 | - 5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) ) & |
---|
398 | + ( sk(k,j,i+2) - sk(k,j,i-3) ) & |
---|
399 | ) * adv_sca_5 |
---|
400 | |
---|
401 | ENDDO |
---|
402 | |
---|
403 | ENDIF |
---|
404 | |
---|
405 | flux_t(0) = 0.0 |
---|
406 | diss_t(0) = 0.0 |
---|
407 | flux_d = 0.0 |
---|
408 | diss_d = 0.0 |
---|
409 | ! |
---|
410 | !-- Now compute the fluxes and tendency terms for the horizontal and |
---|
411 | !-- vertical parts up to the top of the highest topography. |
---|
412 | DO k = nzb+1, nzb_max |
---|
413 | ! |
---|
414 | !-- Note: It is faster to conduct all multiplications explicitly, e.g. |
---|
415 | !-- * adv_sca_5 ... than to determine a factor and multiplicate the |
---|
416 | !-- flux at the end. Indeed, per loop the first requires 6 |
---|
417 | !-- multiplications more, but the the second requires 9 more internal |
---|
418 | !-- calls of IBITS(). |
---|
419 | u_comp = u(k,j,i+1) - u_gtrans |
---|
420 | flux_r(k) = u_comp * ( & |
---|
421 | ( 37.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 & |
---|
422 | + 7.0 * IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3 & |
---|
423 | + IBITS(wall_flags_0(k,j,i),0,1) * adv_sca_1 & |
---|
424 | ) * & |
---|
425 | ( sk(k,j,i+1) + sk(k,j,i) ) & |
---|
426 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 & |
---|
427 | + IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3 & |
---|
428 | ) * & |
---|
429 | ( sk(k,j,i+2) + sk(k,j,i-1) ) & |
---|
430 | + ( IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 & |
---|
431 | ) * ( sk(k,j,i+3) + sk(k,j,i-2) ) & |
---|
432 | ) |
---|
433 | |
---|
434 | diss_r(k) = -ABS( u_comp ) * ( & |
---|
435 | ( 10.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 & |
---|
436 | + 3.0 * IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3 & |
---|
437 | + IBITS(wall_flags_0(k,j,i),0,1) * adv_sca_1 & |
---|
438 | ) * & |
---|
439 | ( sk(k,j,i+1) - sk(k,j,i) ) & |
---|
440 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 & |
---|
441 | + IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3 & |
---|
442 | ) * & |
---|
443 | ( sk(k,j,i+2) - sk(k,j,i-1) ) & |
---|
444 | + ( IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 & |
---|
445 | ) * & |
---|
446 | ( sk(k,j,i+3) - sk(k,j,i-2) ) & |
---|
447 | ) |
---|
448 | |
---|
449 | v_comp = v(k,j+1,i) - v_gtrans |
---|
450 | flux_n(k) = v_comp * ( & |
---|
451 | ( 37.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 & |
---|
452 | + 7.0 * IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3 & |
---|
453 | + IBITS(wall_flags_0(k,j,i),3,1) * adv_sca_1 & |
---|
454 | ) * & |
---|
455 | ( sk(k,j+1,i) + sk(k,j,i) ) & |
---|
456 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 & |
---|
457 | + IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3 & |
---|
458 | ) * & |
---|
459 | ( sk(k,j+2,i) + sk(k,j-1,i) ) & |
---|
460 | + ( IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 & |
---|
461 | ) * & |
---|
462 | ( sk(k,j+3,i) + sk(k,j-2,i) ) & |
---|
463 | ) |
---|
464 | |
---|
465 | diss_n(k) = -ABS( v_comp ) * ( & |
---|
466 | ( 10.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 & |
---|
467 | + 3.0 * IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3 & |
---|
468 | + IBITS(wall_flags_0(k,j,i),3,1) * adv_sca_1 & |
---|
469 | ) * & |
---|
470 | ( sk(k,j+1,i) - sk(k,j,i) ) & |
---|
471 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 & |
---|
472 | + IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3 & |
---|
473 | ) * & |
---|
474 | ( sk(k,j+2,i) - sk(k,j-1,i) ) & |
---|
475 | + ( IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 & |
---|
476 | ) * & |
---|
477 | ( sk(k,j+3,i) - sk(k,j-2,i) ) & |
---|
478 | ) |
---|
479 | ! |
---|
480 | !-- k index has to be modified near bottom and top, else array |
---|
481 | !-- subscripts will be exceeded. |
---|
482 | k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),8,1) |
---|
483 | k_pp = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),6,1) ) |
---|
484 | k_mm = k - 2 * IBITS(wall_flags_0(k,j,i),8,1) |
---|
485 | |
---|
486 | |
---|
487 | flux_t(k) = w(k,j,i) * ( & |
---|
488 | ( 37.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 & |
---|
489 | + 7.0 * IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3 & |
---|
490 | + IBITS(wall_flags_0(k,j,i),6,1) * adv_sca_1 & |
---|
491 | ) * & |
---|
492 | ( sk(k+1,j,i) + sk(k,j,i) ) & |
---|
493 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 & |
---|
494 | + IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3 & |
---|
495 | ) * & |
---|
496 | ( sk(k_pp,j,i) + sk(k-1,j,i) ) & |
---|
497 | + ( IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 & |
---|
498 | ) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) & |
---|
499 | ) |
---|
500 | |
---|
501 | diss_t(k) = -ABS( w(k,j,i) ) * ( & |
---|
502 | ( 10.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 & |
---|
503 | + 3.0 * IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3 & |
---|
504 | + IBITS(wall_flags_0(k,j,i),6,1) * adv_sca_1 & |
---|
505 | ) * & |
---|
506 | ( sk(k+1,j,i) - sk(k,j,i) ) & |
---|
507 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 & |
---|
508 | + IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3 & |
---|
509 | ) * & |
---|
510 | ( sk(k_pp,j,i) - sk(k-1,j,i) ) & |
---|
511 | + ( IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 & |
---|
512 | ) * & |
---|
513 | ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & |
---|
514 | ) |
---|
515 | ! |
---|
516 | !-- Calculate the divergence of the velocity field. A respective |
---|
517 | !-- correction is needed to overcome numerical instabilities caused |
---|
518 | !-- by a not sufficient reduction of divergences near topography. |
---|
519 | div = ( u(k,j,i+1) - u(k,j,i) ) * ddx & |
---|
520 | + ( v(k,j+1,i) - v(k,j,i) ) * ddy & |
---|
521 | + ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) |
---|
522 | |
---|
523 | tend(k,j,i) = tend(k,j,i) - ( & |
---|
524 | ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - & |
---|
525 | swap_diss_x_local(k,j,tn) ) * ddx & |
---|
526 | + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - & |
---|
527 | swap_diss_y_local(k,tn) ) * ddy & |
---|
528 | + ( flux_t(k) + diss_t(k) - flux_d - diss_d & |
---|
529 | ) * ddzw(k) & |
---|
530 | ) + sk(k,j,i) * div |
---|
531 | |
---|
532 | swap_flux_y_local(k,tn) = flux_n(k) |
---|
533 | swap_diss_y_local(k,tn) = diss_n(k) |
---|
534 | swap_flux_x_local(k,j,tn) = flux_r(k) |
---|
535 | swap_diss_x_local(k,j,tn) = diss_r(k) |
---|
536 | flux_d = flux_t(k) |
---|
537 | diss_d = diss_t(k) |
---|
538 | |
---|
539 | ENDDO |
---|
540 | ! |
---|
541 | !-- Now compute the fluxes and tendency terms for the horizontal and |
---|
542 | !-- vertical parts above the top of the highest topography. No degradation |
---|
543 | !-- for the horizontal parts, but for the vertical it is stell needed. |
---|
544 | DO k = nzb_max+1, nzt |
---|
545 | |
---|
546 | u_comp = u(k,j,i+1) - u_gtrans |
---|
547 | flux_r(k) = u_comp * ( & |
---|
548 | 37.0 * ( sk(k,j,i+1) + sk(k,j,i) ) & |
---|
549 | - 8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) ) & |
---|
550 | + ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5 |
---|
551 | diss_r(k) = -ABS( u_comp ) * ( & |
---|
552 | 10.0 * ( sk(k,j,i+1) - sk(k,j,i) ) & |
---|
553 | - 5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) ) & |
---|
554 | + ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5 |
---|
555 | |
---|
556 | v_comp = v(k,j+1,i) - v_gtrans |
---|
557 | flux_n(k) = v_comp * ( & |
---|
558 | 37.0 * ( sk(k,j+1,i) + sk(k,j,i) ) & |
---|
559 | - 8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) ) & |
---|
560 | + ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5 |
---|
561 | diss_n(k) = -ABS( v_comp ) * ( & |
---|
562 | 10.0 * ( sk(k,j+1,i) - sk(k,j,i) ) & |
---|
563 | - 5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) ) & |
---|
564 | + ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5 |
---|
565 | ! |
---|
566 | !-- k index has to be modified near bottom and top, else array |
---|
567 | !-- subscripts will be exceeded. |
---|
568 | k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),8,1) |
---|
569 | k_pp = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),6,1) ) |
---|
570 | k_mm = k - 2 * IBITS(wall_flags_0(k,j,i),8,1) |
---|
571 | |
---|
572 | flux_t(k) = w(k,j,i) * ( & |
---|
573 | ( 37.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 & |
---|
574 | + 7.0 * IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3 & |
---|
575 | + IBITS(wall_flags_0(k,j,i),6,1) * adv_sca_1 & |
---|
576 | ) * & |
---|
577 | ( sk(k+1,j,i) + sk(k,j,i) ) & |
---|
578 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 & |
---|
579 | + IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3 & |
---|
580 | ) * & |
---|
581 | ( sk(k_pp,j,i) + sk(k-1,j,i) ) & |
---|
582 | + ( IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 & |
---|
583 | ) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) & |
---|
584 | ) |
---|
585 | |
---|
586 | diss_t(k) = -ABS( w(k,j,i) ) * ( & |
---|
587 | ( 10.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 & |
---|
588 | + 3.0 * IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3 & |
---|
589 | + IBITS(wall_flags_0(k,j,i),6,1) * adv_sca_1 & |
---|
590 | ) * & |
---|
591 | ( sk(k+1,j,i) - sk(k,j,i) ) & |
---|
592 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 & |
---|
593 | + IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3 & |
---|
594 | ) * & |
---|
595 | ( sk(k_pp,j,i) - sk(k-1,j,i) ) & |
---|
596 | + ( IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 & |
---|
597 | ) * & |
---|
598 | ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & |
---|
599 | ) |
---|
600 | ! |
---|
601 | !-- Calculate the divergence of the velocity field. A respective |
---|
602 | !-- correction is needed to overcome numerical instabilities introduced |
---|
603 | !-- by a not sufficient reduction of divergences near topography. |
---|
604 | div = ( u(k,j,i+1) - u(k,j,i) ) * ddx & |
---|
605 | + ( v(k,j+1,i) - v(k,j,i) ) * ddy & |
---|
606 | + ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) |
---|
607 | |
---|
608 | tend(k,j,i) = tend(k,j,i) - ( & |
---|
609 | ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - & |
---|
610 | swap_diss_x_local(k,j,tn) ) * ddx & |
---|
611 | + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - & |
---|
612 | swap_diss_y_local(k,tn) ) * ddy & |
---|
613 | + ( flux_t(k) + diss_t(k) - flux_d - diss_d & |
---|
614 | ) * ddzw(k) & |
---|
615 | ) + sk(k,j,i) * div |
---|
616 | |
---|
617 | swap_flux_y_local(k,tn) = flux_n(k) |
---|
618 | swap_diss_y_local(k,tn) = diss_n(k) |
---|
619 | swap_flux_x_local(k,j,tn) = flux_r(k) |
---|
620 | swap_diss_x_local(k,j,tn) = diss_r(k) |
---|
621 | flux_d = flux_t(k) |
---|
622 | diss_d = diss_t(k) |
---|
623 | |
---|
624 | ENDDO |
---|
625 | |
---|
626 | ! |
---|
627 | !-- Evaluation of statistics |
---|
628 | SELECT CASE ( sk_char ) |
---|
629 | |
---|
630 | CASE ( 'pt' ) |
---|
631 | |
---|
632 | DO k = nzb, nzt |
---|
633 | sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) + & |
---|
634 | ( flux_t(k) + diss_t(k) ) & |
---|
635 | * weight_substep(intermediate_timestep_count) |
---|
636 | ENDDO |
---|
637 | |
---|
638 | CASE ( 'sa' ) |
---|
639 | |
---|
640 | DO k = nzb, nzt |
---|
641 | sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) + & |
---|
642 | ( flux_t(k) + diss_t(k) ) & |
---|
643 | * weight_substep(intermediate_timestep_count) |
---|
644 | ENDDO |
---|
645 | |
---|
646 | CASE ( 'q' ) |
---|
647 | |
---|
648 | DO k = nzb, nzt |
---|
649 | sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) + & |
---|
650 | ( flux_t(k) + diss_t(k) ) & |
---|
651 | * weight_substep(intermediate_timestep_count) |
---|
652 | ENDDO |
---|
653 | |
---|
654 | END SELECT |
---|
655 | |
---|
656 | END SUBROUTINE advec_s_ws_ij |
---|
657 | |
---|
658 | |
---|
659 | |
---|
660 | |
---|
661 | !------------------------------------------------------------------------------! |
---|
662 | ! Advection of u-component - Call for grid point i,j |
---|
663 | !------------------------------------------------------------------------------! |
---|
664 | SUBROUTINE advec_u_ws_ij( i, j, i_omp, tn ) |
---|
665 | |
---|
666 | USE arrays_3d |
---|
667 | USE constants |
---|
668 | USE control_parameters |
---|
669 | USE grid_variables |
---|
670 | USE indices |
---|
671 | USE statistics |
---|
672 | |
---|
673 | IMPLICIT NONE |
---|
674 | |
---|
675 | INTEGER :: i, i_omp, j, k, k_mm, k_pp, k_ppp, tn |
---|
676 | REAL :: diss_d, div, flux_d, gu, gv, u_comp_l, v_comp, w_comp |
---|
677 | REAL, DIMENSION(nzb:nzt+1) :: diss_n, diss_r, diss_t, flux_n, flux_r, & |
---|
678 | flux_t, u_comp |
---|
679 | |
---|
680 | gu = 2.0 * u_gtrans |
---|
681 | gv = 2.0 * v_gtrans |
---|
682 | ! |
---|
683 | !-- Compute southside fluxes for the respective boundary of PE |
---|
684 | IF ( j == nys ) THEN |
---|
685 | |
---|
686 | DO k = nzb+1, nzb_max |
---|
687 | |
---|
688 | v_comp = v(k,j,i) + v(k,j,i-1) - gv |
---|
689 | flux_s_u(k,tn) = v_comp * ( & |
---|
690 | ( 37.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 & |
---|
691 | + 7.0 * IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 & |
---|
692 | + IBITS(wall_flags_0(k,j,i),12,1) * adv_mom_1 & |
---|
693 | ) * & |
---|
694 | ( u(k,j,i) + u(k,j-1,i) ) & |
---|
695 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 & |
---|
696 | + IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 & |
---|
697 | ) * & |
---|
698 | ( u(k,j+1,i) + u(k,j-2,i) ) & |
---|
699 | + ( IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 & |
---|
700 | ) * & |
---|
701 | ( u(k,j+2,i) + u(k,j-3,i) ) & |
---|
702 | ) |
---|
703 | |
---|
704 | diss_s_u(k,tn) = - ABS ( v_comp ) * ( & |
---|
705 | ( 10.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 & |
---|
706 | + 3.0 * IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 & |
---|
707 | + IBITS(wall_flags_0(k,j,i),12,1) * adv_mom_1 & |
---|
708 | ) * & |
---|
709 | ( u(k,j,i) - u(k,j-1,i) ) & |
---|
710 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 & |
---|
711 | + IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 & |
---|
712 | ) * & |
---|
713 | ( u(k,j+1,i) - u(k,j-2,i) ) & |
---|
714 | + ( IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 & |
---|
715 | ) * & |
---|
716 | ( u(k,j+2,i) - u(k,j-3,i) ) & |
---|
717 | ) |
---|
718 | |
---|
719 | ENDDO |
---|
720 | |
---|
721 | DO k = nzb_max+1, nzt |
---|
722 | |
---|
723 | v_comp = v(k,j,i) + v(k,j,i-1) - gv |
---|
724 | flux_s_u(k,tn) = v_comp * ( & |
---|
725 | 37.0 * ( u(k,j,i) + u(k,j-1,i) ) & |
---|
726 | - 8.0 * ( u(k,j+1,i) + u(k,j-2,i) ) & |
---|
727 | + ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5 |
---|
728 | diss_s_u(k,tn) = - ABS(v_comp) * ( & |
---|
729 | 10.0 * ( u(k,j,i) - u(k,j-1,i) ) & |
---|
730 | - 5.0 * ( u(k,j+1,i) - u(k,j-2,i) ) & |
---|
731 | + ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5 |
---|
732 | |
---|
733 | ENDDO |
---|
734 | |
---|
735 | ENDIF |
---|
736 | ! |
---|
737 | !-- Compute leftside fluxes for the respective boundary of PE |
---|
738 | IF ( i == i_omp ) THEN |
---|
739 | |
---|
740 | DO k = nzb+1, nzb_max |
---|
741 | |
---|
742 | u_comp_l = u(k,j,i) + u(k,j,i-1) - gu |
---|
743 | flux_l_u(k,j,tn) = u_comp_l * ( & |
---|
744 | ( 37.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 & |
---|
745 | + 7.0 * IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 & |
---|
746 | + IBITS(wall_flags_0(k,j,i),9,1) * adv_mom_1 & |
---|
747 | ) * & |
---|
748 | ( u(k,j,i) + u(k,j,i-1) ) & |
---|
749 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 & |
---|
750 | + IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 & |
---|
751 | ) * & |
---|
752 | ( u(k,j,i+1) + u(k,j,i-2) ) & |
---|
753 | + ( IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 & |
---|
754 | ) * & |
---|
755 | ( u(k,j,i+2) + u(k,j,i-3) ) & |
---|
756 | ) |
---|
757 | |
---|
758 | diss_l_u(k,j,tn) = - ABS( u_comp_l ) * ( & |
---|
759 | ( 10.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 & |
---|
760 | + 3.0 * IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 & |
---|
761 | + IBITS(wall_flags_0(k,j,i),9,1) * adv_mom_1 & |
---|
762 | ) * & |
---|
763 | ( u(k,j,i) - u(k,j,i-1) ) & |
---|
764 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 & |
---|
765 | + IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 & |
---|
766 | ) * & |
---|
767 | ( u(k,j,i+1) - u(k,j,i-2) ) & |
---|
768 | + ( IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 & |
---|
769 | ) * & |
---|
770 | ( u(k,j,i+2) - u(k,j,i-3) ) & |
---|
771 | ) |
---|
772 | |
---|
773 | ENDDO |
---|
774 | |
---|
775 | DO k = nzb_max+1, nzt |
---|
776 | |
---|
777 | u_comp_l = u(k,j,i) + u(k,j,i-1) - gu |
---|
778 | flux_l_u(k,j,tn) = u_comp_l * ( & |
---|
779 | 37.0 * ( u(k,j,i) + u(k,j,i-1) ) & |
---|
780 | - 8.0 * ( u(k,j,i+1) + u(k,j,i-2) ) & |
---|
781 | + ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5 |
---|
782 | diss_l_u(k,j,tn) = - ABS(u_comp_l) * ( & |
---|
783 | 10.0 * ( u(k,j,i) - u(k,j,i-1) ) & |
---|
784 | - 5.0 * ( u(k,j,i+1) - u(k,j,i-2) ) & |
---|
785 | + ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5 |
---|
786 | |
---|
787 | ENDDO |
---|
788 | |
---|
789 | ENDIF |
---|
790 | |
---|
791 | flux_t(0) = 0.0 |
---|
792 | diss_t(0) = 0.0 |
---|
793 | flux_d = 0.0 |
---|
794 | diss_d = 0.0 |
---|
795 | ! |
---|
796 | !-- Now compute the fluxes tendency terms for the horizontal and |
---|
797 | !-- vertical parts. |
---|
798 | DO k = nzb+1, nzb_max |
---|
799 | |
---|
800 | u_comp(k) = u(k,j,i+1) + u(k,j,i) |
---|
801 | flux_r(k) = ( u_comp(k) - gu ) * ( & |
---|
802 | ( 37.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 & |
---|
803 | + 7.0 * IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 & |
---|
804 | + IBITS(wall_flags_0(k,j,i),9,1) * adv_mom_1 & |
---|
805 | ) * & |
---|
806 | ( u(k,j,i+1) + u(k,j,i) ) & |
---|
807 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 & |
---|
808 | + IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 & |
---|
809 | ) * & |
---|
810 | ( u(k,j,i+2) + u(k,j,i-1) ) & |
---|
811 | + ( IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 & |
---|
812 | ) * & |
---|
813 | ( u(k,j,i+3) + u(k,j,i-2) ) & |
---|
814 | ) |
---|
815 | |
---|
816 | diss_r(k) = - ABS( u_comp(k) - gu ) * ( & |
---|
817 | ( 10.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 & |
---|
818 | + 3.0 * IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 & |
---|
819 | + IBITS(wall_flags_0(k,j,i),9,1) * adv_mom_1 & |
---|
820 | ) * & |
---|
821 | ( u(k,j,i+1) - u(k,j,i) ) & |
---|
822 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 & |
---|
823 | + IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 & |
---|
824 | ) * & |
---|
825 | ( u(k,j,i+2) - u(k,j,i-1) ) & |
---|
826 | + ( IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 & |
---|
827 | ) * & |
---|
828 | ( u(k,j,i+3) - u(k,j,i-2) ) & |
---|
829 | ) |
---|
830 | |
---|
831 | v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv |
---|
832 | flux_n(k) = v_comp * ( & |
---|
833 | ( 37.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 & |
---|
834 | + 7.0 * IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 & |
---|
835 | + IBITS(wall_flags_0(k,j,i),12,1) * adv_mom_1 & |
---|
836 | ) * & |
---|
837 | ( u(k,j+1,i) + u(k,j,i) ) & |
---|
838 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 & |
---|
839 | + IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 & |
---|
840 | ) * & |
---|
841 | ( u(k,j+2,i) + u(k,j-1,i) ) & |
---|
842 | + ( IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 & |
---|
843 | ) * & |
---|
844 | ( u(k,j+3,i) + u(k,j-2,i) ) & |
---|
845 | ) |
---|
846 | |
---|
847 | diss_n(k) = - ABS ( v_comp ) * ( & |
---|
848 | ( 10.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 & |
---|
849 | + 3.0 * IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 & |
---|
850 | + IBITS(wall_flags_0(k,j,i),12,1) * adv_mom_1 & |
---|
851 | ) * & |
---|
852 | ( u(k,j+1,i) - u(k,j,i) ) & |
---|
853 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 & |
---|
854 | + IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 & |
---|
855 | ) * & |
---|
856 | ( u(k,j+2,i) - u(k,j-1,i) ) & |
---|
857 | + ( IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 & |
---|
858 | ) * & |
---|
859 | ( u(k,j+3,i) - u(k,j-2,i) ) & |
---|
860 | ) |
---|
861 | ! |
---|
862 | !-- k index has to be modified near bottom and top, else array |
---|
863 | !-- subscripts will be exceeded. |
---|
864 | k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),17,1) |
---|
865 | k_pp = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),15,1) ) |
---|
866 | k_mm = k - 2 * IBITS(wall_flags_0(k,j,i),17,1) |
---|
867 | |
---|
868 | w_comp = w(k,j,i) + w(k,j,i-1) |
---|
869 | flux_t(k) = w_comp * ( & |
---|
870 | ( 37.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 & |
---|
871 | + 7.0 * IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3 & |
---|
872 | + IBITS(wall_flags_0(k,j,i),15,1) * adv_mom_1 & |
---|
873 | ) * & |
---|
874 | ( u(k+1,j,i) + u(k,j,i) ) & |
---|
875 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 & |
---|
876 | + IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3 & |
---|
877 | ) * & |
---|
878 | ( u(k_pp,j,i) + u(k-1,j,i) ) & |
---|
879 | + ( IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 & |
---|
880 | ) * & |
---|
881 | ( u(k_ppp,j,i) + u(k_mm,j,i) ) & |
---|
882 | ) |
---|
883 | |
---|
884 | diss_t(k) = - ABS( w_comp ) * ( & |
---|
885 | ( 10.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 & |
---|
886 | + 3.0 * IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3 & |
---|
887 | + IBITS(wall_flags_0(k,j,i),15,1) * adv_mom_1 & |
---|
888 | ) * & |
---|
889 | ( u(k+1,j,i) - u(k,j,i) ) & |
---|
890 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 & |
---|
891 | + IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3 & |
---|
892 | ) * & |
---|
893 | ( u(k_pp,j,i) - u(k-1,j,i) ) & |
---|
894 | + ( IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 & |
---|
895 | ) * & |
---|
896 | ( u(k_ppp,j,i) - u(k_mm,j,i) ) & |
---|
897 | ) |
---|
898 | ! |
---|
899 | !-- Calculate the divergence of the velocity field. A respective |
---|
900 | !-- correction is needed to overcome numerical instabilities introduced |
---|
901 | !-- by a not sufficient reduction of divergences near topography. |
---|
902 | div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx & |
---|
903 | + ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy & |
---|
904 | + ( w_comp - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) * ddzw(k) & |
---|
905 | ) * 0.5 |
---|
906 | |
---|
907 | tend(k,j,i) = tend(k,j,i) - ( & |
---|
908 | ( flux_r(k) + diss_r(k) & |
---|
909 | - flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx & |
---|
910 | + ( flux_n(k) + diss_n(k) & |
---|
911 | - flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy & |
---|
912 | + ( flux_t(k) + diss_t(k) & |
---|
913 | - flux_d - diss_d ) * ddzw(k) & |
---|
914 | ) + div * u(k,j,i) |
---|
915 | |
---|
916 | flux_l_u(k,j,tn) = flux_r(k) |
---|
917 | diss_l_u(k,j,tn) = diss_r(k) |
---|
918 | flux_s_u(k,tn) = flux_n(k) |
---|
919 | diss_s_u(k,tn) = diss_n(k) |
---|
920 | flux_d = flux_t(k) |
---|
921 | diss_d = diss_t(k) |
---|
922 | ! |
---|
923 | !-- Statistical Evaluation of u'u'. The factor has to be applied for |
---|
924 | !-- right evaluation when gallilei_trans = .T. . |
---|
925 | sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) & |
---|
926 | + ( flux_r(k) * & |
---|
927 | ( u_comp(k) - 2.0 * hom(k,1,1,0) ) & |
---|
928 | / ( u_comp(k) - gu + 1.0E-20 ) & |
---|
929 | + diss_r(k) * & |
---|
930 | ABS( u_comp(k) - 2.0 * hom(k,1,1,0) ) & |
---|
931 | / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) ) & |
---|
932 | * weight_substep(intermediate_timestep_count) |
---|
933 | ! |
---|
934 | !-- Statistical Evaluation of w'u'. |
---|
935 | sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) & |
---|
936 | + ( flux_t(k) + diss_t(k) ) & |
---|
937 | * weight_substep(intermediate_timestep_count) |
---|
938 | ENDDO |
---|
939 | |
---|
940 | DO k = nzb_max+1, nzt |
---|
941 | |
---|
942 | u_comp(k) = u(k,j,i+1) + u(k,j,i) |
---|
943 | flux_r(k) = ( u_comp(k) - gu ) * ( & |
---|
944 | 37.0 * ( u(k,j,i+1) + u(k,j,i) ) & |
---|
945 | - 8.0 * ( u(k,j,i+2) + u(k,j,i-1) ) & |
---|
946 | + ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5 |
---|
947 | diss_r(k) = - ABS( u_comp(k) - gu ) * ( & |
---|
948 | 10.0 * ( u(k,j,i+1) - u(k,j,i) ) & |
---|
949 | - 5.0 * ( u(k,j,i+2) - u(k,j,i-1) ) & |
---|
950 | + ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5 |
---|
951 | |
---|
952 | v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv |
---|
953 | flux_n(k) = v_comp * ( & |
---|
954 | 37.0 * ( u(k,j+1,i) + u(k,j,i) ) & |
---|
955 | - 8.0 * ( u(k,j+2,i) + u(k,j-1,i) ) & |
---|
956 | + ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5 |
---|
957 | diss_n(k) = - ABS( v_comp ) * ( & |
---|
958 | 10.0 * ( u(k,j+1,i) - u(k,j,i) ) & |
---|
959 | - 5.0 * ( u(k,j+2,i) - u(k,j-1,i) ) & |
---|
960 | + ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5 |
---|
961 | ! |
---|
962 | !-- k index has to be modified near bottom and top, else array |
---|
963 | !-- subscripts will be exceeded. |
---|
964 | k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),17,1) |
---|
965 | k_pp = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),15,1) ) |
---|
966 | k_mm = k - 2 * IBITS(wall_flags_0(k,j,i),17,1) |
---|
967 | |
---|
968 | w_comp = w(k,j,i) + w(k,j,i-1) |
---|
969 | flux_t(k) = w_comp * ( & |
---|
970 | ( 37.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 & |
---|
971 | + 7.0 * IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3 & |
---|
972 | + IBITS(wall_flags_0(k,j,i),15,1) * adv_mom_1 & |
---|
973 | ) * & |
---|
974 | ( u(k+1,j,i) + u(k,j,i) ) & |
---|
975 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 & |
---|
976 | + IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3 & |
---|
977 | ) * & |
---|
978 | ( u(k_pp,j,i) + u(k-1,j,i) ) & |
---|
979 | + ( IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 & |
---|
980 | ) * & |
---|
981 | ( u(k_ppp,j,i) + u(k_mm,j,i) ) & |
---|
982 | ) |
---|
983 | |
---|
984 | diss_t(k) = - ABS( w_comp ) * ( & |
---|
985 | ( 10.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 & |
---|
986 | + 3.0 * IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3 & |
---|
987 | + IBITS(wall_flags_0(k,j,i),15,1) * adv_mom_1 & |
---|
988 | ) * & |
---|
989 | ( u(k+1,j,i) - u(k,j,i) ) & |
---|
990 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 & |
---|
991 | + IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3 & |
---|
992 | ) * & |
---|
993 | ( u(k_pp,j,i) - u(k-1,j,i) ) & |
---|
994 | + ( IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 & |
---|
995 | ) * & |
---|
996 | ( u(k_ppp,j,i) - u(k_mm,j,i) ) & |
---|
997 | ) |
---|
998 | ! |
---|
999 | !-- Calculate the divergence of the velocity field. A respective |
---|
1000 | !-- correction is needed to overcome numerical instabilities introduced |
---|
1001 | !-- by a not sufficient reduction of divergences near topography. |
---|
1002 | div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx & |
---|
1003 | + ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy & |
---|
1004 | + ( w_comp - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) * ddzw(k) & |
---|
1005 | ) * 0.5 |
---|
1006 | |
---|
1007 | tend(k,j,i) = tend(k,j,i) - ( & |
---|
1008 | ( flux_r(k) + diss_r(k) & |
---|
1009 | - flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx & |
---|
1010 | + ( flux_n(k) + diss_n(k) & |
---|
1011 | - flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy & |
---|
1012 | + ( flux_t(k) + diss_t(k) & |
---|
1013 | - flux_d - diss_d ) * ddzw(k) & |
---|
1014 | ) + div * u(k,j,i) |
---|
1015 | |
---|
1016 | flux_l_u(k,j,tn) = flux_r(k) |
---|
1017 | diss_l_u(k,j,tn) = diss_r(k) |
---|
1018 | flux_s_u(k,tn) = flux_n(k) |
---|
1019 | diss_s_u(k,tn) = diss_n(k) |
---|
1020 | flux_d = flux_t(k) |
---|
1021 | diss_d = diss_t(k) |
---|
1022 | ! |
---|
1023 | !-- Statistical Evaluation of u'u'. The factor has to be applied for |
---|
1024 | !-- right evaluation when gallilei_trans = .T. . |
---|
1025 | sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) & |
---|
1026 | + ( flux_r(k) * & |
---|
1027 | ( u_comp(k) - 2.0 * hom(k,1,1,0) ) & |
---|
1028 | / ( u_comp(k) - gu + 1.0E-20 ) & |
---|
1029 | + diss_r(k) * & |
---|
1030 | ABS( u_comp(k) - 2.0 * hom(k,1,1,0) ) & |
---|
1031 | / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) ) & |
---|
1032 | * weight_substep(intermediate_timestep_count) |
---|
1033 | ! |
---|
1034 | !-- Statistical Evaluation of w'u'. |
---|
1035 | sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) & |
---|
1036 | + ( flux_t(k) + diss_t(k) ) & |
---|
1037 | * weight_substep(intermediate_timestep_count) |
---|
1038 | ENDDO |
---|
1039 | |
---|
1040 | sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn) |
---|
1041 | |
---|
1042 | |
---|
1043 | |
---|
1044 | END SUBROUTINE advec_u_ws_ij |
---|
1045 | |
---|
1046 | |
---|
1047 | |
---|
1048 | !-----------------------------------------------------------------------------! |
---|
1049 | ! Advection of v-component - Call for grid point i,j |
---|
1050 | !-----------------------------------------------------------------------------! |
---|
1051 | SUBROUTINE advec_v_ws_ij( i, j, i_omp, tn ) |
---|
1052 | |
---|
1053 | USE arrays_3d |
---|
1054 | USE constants |
---|
1055 | USE control_parameters |
---|
1056 | USE grid_variables |
---|
1057 | USE indices |
---|
1058 | USE statistics |
---|
1059 | |
---|
1060 | IMPLICIT NONE |
---|
1061 | |
---|
1062 | INTEGER :: i, i_omp, j, k, k_mm, k_pp, k_ppp, tn |
---|
1063 | REAL :: diss_d, div, flux_d, gu, gv, u_comp, v_comp_l, w_comp |
---|
1064 | REAL, DIMENSION(nzb:nzt+1) :: diss_n, diss_r, diss_t, flux_n, flux_r, & |
---|
1065 | flux_t, v_comp |
---|
1066 | |
---|
1067 | gu = 2.0 * u_gtrans |
---|
1068 | gv = 2.0 * v_gtrans |
---|
1069 | |
---|
1070 | ! |
---|
1071 | !-- Compute leftside fluxes for the respective boundary. |
---|
1072 | IF ( i == i_omp ) THEN |
---|
1073 | |
---|
1074 | DO k = nzb+1, nzb_max |
---|
1075 | |
---|
1076 | u_comp = u(k,j-1,i) + u(k,j,i) - gu |
---|
1077 | flux_l_v(k,j,tn) = u_comp * ( & |
---|
1078 | ( 37.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 & |
---|
1079 | + 7.0 * IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3 & |
---|
1080 | + IBITS(wall_flags_0(k,j,i),18,1) * adv_mom_1 & |
---|
1081 | ) * & |
---|
1082 | ( v(k,j,i) + v(k,j,i-1) ) & |
---|
1083 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 & |
---|
1084 | + IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3 & |
---|
1085 | ) * & |
---|
1086 | ( v(k,j,i+1) + v(k,j,i-2) ) & |
---|
1087 | + ( IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 & |
---|
1088 | ) * & |
---|
1089 | ( v(k,j,i+2) + v(k,j,i-3) ) & |
---|
1090 | ) |
---|
1091 | |
---|
1092 | diss_l_v(k,j,tn) = - ABS( u_comp ) * ( & |
---|
1093 | ( 10.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 & |
---|
1094 | + 3.0 * IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3 & |
---|
1095 | + IBITS(wall_flags_0(k,j,i),18,1) * adv_mom_1 & |
---|
1096 | ) * & |
---|
1097 | ( v(k,j,i) - v(k,j,i-1) ) & |
---|
1098 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 & |
---|
1099 | + IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3 & |
---|
1100 | ) * & |
---|
1101 | ( v(k,j,i+1) - v(k,j,i-2) ) & |
---|
1102 | + ( IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 & |
---|
1103 | ) * & |
---|
1104 | ( v(k,j,i+2) - v(k,j,i-3) ) & |
---|
1105 | ) |
---|
1106 | |
---|
1107 | ENDDO |
---|
1108 | |
---|
1109 | DO k = nzb_max+1, nzt |
---|
1110 | |
---|
1111 | u_comp = u(k,j-1,i) + u(k,j,i) - gu |
---|
1112 | flux_l_v(k,j,tn) = u_comp * ( & |
---|
1113 | 37.0 * ( v(k,j,i) + v(k,j,i-1) ) & |
---|
1114 | - 8.0 * ( v(k,j,i+1) + v(k,j,i-2) ) & |
---|
1115 | + ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5 |
---|
1116 | diss_l_v(k,j,tn) = - ABS( u_comp ) * ( & |
---|
1117 | 10.0 * ( v(k,j,i) - v(k,j,i-1) ) & |
---|
1118 | - 5.0 * ( v(k,j,i+1) - v(k,j,i-2) ) & |
---|
1119 | + ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5 |
---|
1120 | |
---|
1121 | ENDDO |
---|
1122 | |
---|
1123 | ENDIF |
---|
1124 | ! |
---|
1125 | !-- Compute southside fluxes for the respective boundary. |
---|
1126 | IF ( j == nysv ) THEN |
---|
1127 | |
---|
1128 | DO k = nzb+1, nzb_max |
---|
1129 | |
---|
1130 | v_comp_l = v(k,j,i) + v(k,j-1,i) - gv |
---|
1131 | flux_s_v(k,tn) = v_comp_l * ( & |
---|
1132 | ( 37.0 * IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5 & |
---|
1133 | + 7.0 * IBITS(wall_flags_0(k,j,i),22,1) * adv_mom_3 & |
---|
1134 | + IBITS(wall_flags_0(k,j,i),21,1) * adv_mom_1 & |
---|
1135 | ) * & |
---|
1136 | ( v(k,j,i) + v(k,j-1,i) ) & |
---|
1137 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5 & |
---|
1138 | + IBITS(wall_flags_0(k,j,i),22,1) * adv_mom_3 & |
---|
1139 | ) * & |
---|
1140 | ( v(k,j+1,i) + v(k,j-2,i) ) & |
---|
1141 | + ( IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5 & |
---|
1142 | ) * & |
---|
1143 | ( v(k,j+2,i) + v(k,j-3,i) ) & |
---|
1144 | ) |
---|
1145 | |
---|
1146 | diss_s_v(k,tn) = - ABS( v_comp_l ) * ( & |
---|
1147 | ( 10.0 * IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5 & |
---|
1148 | + 3.0 * IBITS(wall_flags_0(k,j,i),22,1) * adv_mom_3 & |
---|
1149 | + IBITS(wall_flags_0(k,j,i),21,1) * adv_mom_1 & |
---|
1150 | ) * & |
---|
1151 | ( v(k,j,i) - v(k,j-1,i) ) & |
---|
1152 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5 & |
---|
1153 | + IBITS(wall_flags_0(k,j,i),22,1) * adv_mom_3 & |
---|
1154 | ) * & |
---|
1155 | ( v(k,j+1,i) - v(k,j-2,i) ) & |
---|
1156 | + ( IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5 & |
---|
1157 | ) * & |
---|
1158 | ( v(k,j+2,i) - v(k,j-3,i) ) & |
---|
1159 | ) |
---|
1160 | |
---|
1161 | ENDDO |
---|
1162 | |
---|
1163 | DO k = nzb_max+1, nzt |
---|
1164 | |
---|
1165 | v_comp_l = v(k,j,i) + v(k,j-1,i) - gv |
---|
1166 | flux_s_v(k,tn) = v_comp_l * ( & |
---|
1167 | 37.0 * ( v(k,j,i) + v(k,j-1,i) ) & |
---|
1168 | - 8.0 * ( v(k,j+1,i) + v(k,j-2,i) ) & |
---|
1169 | + ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5 |
---|
1170 | diss_s_v(k,tn) = - ABS( v_comp_l ) * ( & |
---|
1171 | 10.0 * ( v(k,j,i) - v(k,j-1,i) ) & |
---|
1172 | - 5.0 * ( v(k,j+1,i) - v(k,j-2,i) ) & |
---|
1173 | + ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5 |
---|
1174 | |
---|
1175 | ENDDO |
---|
1176 | |
---|
1177 | ENDIF |
---|
1178 | |
---|
1179 | flux_t(0) = 0.0 |
---|
1180 | diss_t(0) = 0.0 |
---|
1181 | flux_d = 0.0 |
---|
1182 | diss_d = 0.0 |
---|
1183 | ! |
---|
1184 | !-- Now compute the fluxes and tendency terms for the horizontal and |
---|
1185 | !-- verical parts. |
---|
1186 | DO k = nzb+1, nzb_max |
---|
1187 | |
---|
1188 | u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu |
---|
1189 | flux_r(k) = u_comp * ( & |
---|
1190 | ( 37.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 & |
---|
1191 | + 7.0 * IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3 & |
---|
1192 | + IBITS(wall_flags_0(k,j,i),18,1) * adv_mom_1 & |
---|
1193 | ) * & |
---|
1194 | ( v(k,j,i+1) + v(k,j,i) ) & |
---|
1195 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 & |
---|
1196 | + IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3 & |
---|
1197 | ) * & |
---|
1198 | ( v(k,j,i+2) + v(k,j,i-1) ) & |
---|
1199 | + ( IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 & |
---|
1200 | ) * & |
---|
1201 | ( v(k,j,i+3) + v(k,j,i-2) ) & |
---|
1202 | ) |
---|
1203 | |
---|
1204 | diss_r(k) = - ABS( u_comp ) * ( & |
---|
1205 | ( 10.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 & |
---|
1206 | + 3.0 * IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3 & |
---|
1207 | + IBITS(wall_flags_0(k,j,i),18,1) * adv_mom_1 & |
---|
1208 | ) * & |
---|
1209 | ( v(k,j,i+1) - v(k,j,i) ) & |
---|
1210 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 & |
---|
1211 | + IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3 & |
---|
1212 | ) * & |
---|
1213 | ( v(k,j,i+2) - v(k,j,i-1) ) & |
---|
1214 | + ( IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 & |
---|
1215 | ) * & |
---|
1216 | ( v(k,j,i+3) - v(k,j,i-2) ) & |
---|
1217 | ) |
---|
1218 | |
---|
1219 | v_comp(k) = v(k,j+1,i) + v(k,j,i) |
---|
1220 | flux_n(k) = ( v_comp(k) - gv ) * ( & |
---|
1221 | ( 37.0 * IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5 & |
---|
1222 | + 7.0 * IBITS(wall_flags_0(k,j,i),22,1) * adv_mom_3 & |
---|
1223 | + IBITS(wall_flags_0(k,j,i),21,1) * adv_mom_1 & |
---|
1224 | ) * & |
---|
1225 | ( v(k,j+1,i) + v(k,j,i) ) & |
---|
1226 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5 & |
---|
1227 | + IBITS(wall_flags_0(k,j,i),22,1) * adv_mom_3 & |
---|
1228 | ) * & |
---|
1229 | ( v(k,j+2,i) + v(k,j-1,i) ) & |
---|
1230 | + ( IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5 & |
---|
1231 | ) * & |
---|
1232 | ( v(k,j+3,i) + v(k,j-2,i) ) & |
---|
1233 | ) |
---|
1234 | |
---|
1235 | diss_n(k) = - ABS( v_comp(k) - gv ) * ( & |
---|
1236 | ( 10.0 * IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5 & |
---|
1237 | + 3.0 * IBITS(wall_flags_0(k,j,i),22,1) * adv_mom_3 & |
---|
1238 | + IBITS(wall_flags_0(k,j,i),21,1) * adv_mom_1 & |
---|
1239 | ) * & |
---|
1240 | ( v(k,j+1,i) - v(k,j,i) ) & |
---|
1241 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5 & |
---|
1242 | + IBITS(wall_flags_0(k,j,i),22,1) * adv_mom_3 & |
---|
1243 | ) * & |
---|
1244 | ( v(k,j+2,i) - v(k,j-1,i) ) & |
---|
1245 | + ( IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5 & |
---|
1246 | ) * & |
---|
1247 | ( v(k,j+3,i) - v(k,j-2,i) ) & |
---|
1248 | ) |
---|
1249 | ! |
---|
1250 | !-- k index has to be modified near bottom and top, else array |
---|
1251 | !-- subscripts will be exceeded. |
---|
1252 | k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),26,1) |
---|
1253 | k_pp = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),24,1) ) |
---|
1254 | k_mm = k - 2 * IBITS(wall_flags_0(k,j,i),26,1) |
---|
1255 | |
---|
1256 | w_comp = w(k,j-1,i) + w(k,j,i) |
---|
1257 | flux_t(k) = w_comp * ( & |
---|
1258 | ( 37.0 * IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5 & |
---|
1259 | + 7.0 * IBITS(wall_flags_0(k,j,i),25,1) * adv_mom_3 & |
---|
1260 | + IBITS(wall_flags_0(k,j,i),24,1) * adv_mom_1 & |
---|
1261 | ) * & |
---|
1262 | ( v(k+1,j,i) + v(k,j,i) ) & |
---|
1263 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5 & |
---|
1264 | + IBITS(wall_flags_0(k,j,i),25,1) * adv_mom_3 & |
---|
1265 | ) * & |
---|
1266 | ( v(k_pp,j,i) + v(k-1,j,i) ) & |
---|
1267 | + ( IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5 & |
---|
1268 | ) * & |
---|
1269 | ( v(k_ppp,j,i) + v(k_mm,j,i) ) & |
---|
1270 | ) |
---|
1271 | |
---|
1272 | diss_t(k) = - ABS( w_comp ) * ( & |
---|
1273 | ( 10.0 * IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5 & |
---|
1274 | + 3.0 * IBITS(wall_flags_0(k,j,i),25,1) * adv_mom_3 & |
---|
1275 | + IBITS(wall_flags_0(k,j,i),24,1) * adv_mom_1 & |
---|
1276 | ) * & |
---|
1277 | ( v(k+1,j,i) - v(k,j,i) ) & |
---|
1278 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5 & |
---|
1279 | + IBITS(wall_flags_0(k,j,i),25,1) * adv_mom_3 & |
---|
1280 | ) * & |
---|
1281 | ( v(k_pp,j,i) - v(k-1,j,i) ) & |
---|
1282 | + ( IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5 & |
---|
1283 | ) * & |
---|
1284 | ( v(k_ppp,j,i) - v(k_mm,j,i) ) & |
---|
1285 | ) |
---|
1286 | ! |
---|
1287 | !-- Calculate the divergence of the velocity field. A respective |
---|
1288 | !-- correction is needed to overcome numerical instabilities introduced |
---|
1289 | !-- by a not sufficient reduction of divergences near topography. |
---|
1290 | div = ( ( u_comp + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx & |
---|
1291 | + ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy & |
---|
1292 | + ( w_comp - ( w(k-1,j-1,i) + w(k-1,j,i) ) ) * ddzw(k) & |
---|
1293 | ) * 0.5 |
---|
1294 | |
---|
1295 | tend(k,j,i) = tend(k,j,i) - ( & |
---|
1296 | ( flux_r(k) + diss_r(k) & |
---|
1297 | - flux_l_v(k,j,tn) - diss_l_v(k,j,tn) ) * ddx & |
---|
1298 | + ( flux_n(k) + diss_n(k) & |
---|
1299 | - flux_s_v(k,tn) - diss_s_v(k,tn) ) * ddy & |
---|
1300 | + ( flux_t(k) + diss_t(k) & |
---|
1301 | - flux_d - diss_d ) * ddzw(k) & |
---|
1302 | ) + v(k,j,i) * div |
---|
1303 | |
---|
1304 | flux_l_v(k,j,tn) = flux_r(k) |
---|
1305 | diss_l_v(k,j,tn) = diss_r(k) |
---|
1306 | flux_s_v(k,tn) = flux_n(k) |
---|
1307 | diss_s_v(k,tn) = diss_n(k) |
---|
1308 | flux_d = flux_t(k) |
---|
1309 | diss_d = diss_t(k) |
---|
1310 | |
---|
1311 | ! |
---|
1312 | !-- Statistical Evaluation of v'v'. The factor has to be applied for |
---|
1313 | !-- right evaluation when gallilei_trans = .T. . |
---|
1314 | sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) & |
---|
1315 | + ( flux_n(k) & |
---|
1316 | * ( v_comp(k) - 2.0 * hom(k,1,2,0) ) & |
---|
1317 | / ( v_comp(k) - gv + 1.0E-20 ) & |
---|
1318 | + diss_n(k) & |
---|
1319 | * ABS( v_comp(k) - 2.0 * hom(k,1,2,0) ) & |
---|
1320 | / ( ABS( v_comp(k) - gv ) +1.0E-20 ) ) & |
---|
1321 | * weight_substep(intermediate_timestep_count) |
---|
1322 | ! |
---|
1323 | !-- Statistical Evaluation of w'v'. |
---|
1324 | sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) & |
---|
1325 | + ( flux_t(k) + diss_t(k) ) & |
---|
1326 | * weight_substep(intermediate_timestep_count) |
---|
1327 | |
---|
1328 | ENDDO |
---|
1329 | |
---|
1330 | DO k = nzb_max+1, nzt |
---|
1331 | |
---|
1332 | u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu |
---|
1333 | flux_r(k) = u_comp * ( & |
---|
1334 | 37.0 * ( v(k,j,i+1) + v(k,j,i) ) & |
---|
1335 | - 8.0 * ( v(k,j,i+2) + v(k,j,i-1) ) & |
---|
1336 | + ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5 |
---|
1337 | |
---|
1338 | diss_r(k) = - ABS( u_comp ) * ( & |
---|
1339 | 10.0 * ( v(k,j,i+1) - v(k,j,i) ) & |
---|
1340 | - 5.0 * ( v(k,j,i+2) - v(k,j,i-1) ) & |
---|
1341 | + ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5 |
---|
1342 | |
---|
1343 | |
---|
1344 | v_comp(k) = v(k,j+1,i) + v(k,j,i) |
---|
1345 | flux_n(k) = ( v_comp(k) - gv ) * ( & |
---|
1346 | 37.0 * ( v(k,j+1,i) + v(k,j,i) ) & |
---|
1347 | - 8.0 * ( v(k,j+2,i) + v(k,j-1,i) ) & |
---|
1348 | + ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5 |
---|
1349 | |
---|
1350 | diss_n(k) = - ABS( v_comp(k) - gv ) * ( & |
---|
1351 | 10.0 * ( v(k,j+1,i) - v(k,j,i) ) & |
---|
1352 | - 5.0 * ( v(k,j+2,i) - v(k,j-1,i) ) & |
---|
1353 | + ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5 |
---|
1354 | ! |
---|
1355 | !-- k index has to be modified near bottom and top, else array |
---|
1356 | !-- subscripts will be exceeded. |
---|
1357 | k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),26,1) |
---|
1358 | k_pp = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),24,1) ) |
---|
1359 | k_mm = k - 2 * IBITS(wall_flags_0(k,j,i),26,1) |
---|
1360 | |
---|
1361 | w_comp = w(k,j-1,i) + w(k,j,i) |
---|
1362 | flux_t(k) = w_comp * ( & |
---|
1363 | ( 37.0 * IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5 & |
---|
1364 | + 7.0 * IBITS(wall_flags_0(k,j,i),25,1) * adv_mom_3 & |
---|
1365 | + IBITS(wall_flags_0(k,j,i),24,1) * adv_mom_1 & |
---|
1366 | ) * & |
---|
1367 | ( v(k+1,j,i) + v(k,j,i) ) & |
---|
1368 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5 & |
---|
1369 | + IBITS(wall_flags_0(k,j,i),25,1) * adv_mom_3 & |
---|
1370 | ) * & |
---|
1371 | ( v(k_pp,j,i) + v(k-1,j,i) ) & |
---|
1372 | + ( IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5 & |
---|
1373 | ) * & |
---|
1374 | ( v(k_ppp,j,i) + v(k_mm,j,i) ) & |
---|
1375 | ) |
---|
1376 | |
---|
1377 | diss_t(k) = - ABS( w_comp ) * ( & |
---|
1378 | ( 10.0 * IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5 & |
---|
1379 | + 3.0 * IBITS(wall_flags_0(k,j,i),25,1) * adv_mom_3 & |
---|
1380 | + IBITS(wall_flags_0(k,j,i),24,1) * adv_mom_1 & |
---|
1381 | ) * & |
---|
1382 | ( v(k+1,j,i) - v(k,j,i) ) & |
---|
1383 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5 & |
---|
1384 | + IBITS(wall_flags_0(k,j,i),25,1) * adv_mom_3 & |
---|
1385 | ) * & |
---|
1386 | ( v(k_pp,j,i) - v(k-1,j,i) ) & |
---|
1387 | + ( IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5 & |
---|
1388 | ) * & |
---|
1389 | ( v(k_ppp,j,i) - v(k_mm,j,i) ) & |
---|
1390 | ) |
---|
1391 | ! |
---|
1392 | !-- Calculate the divergence of the velocity field. A respective |
---|
1393 | !-- correction is needed to overcome numerical instabilities introduced |
---|
1394 | !-- by a not sufficient reduction of divergences near topography. |
---|
1395 | div = ( ( u_comp + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx & |
---|
1396 | + ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy & |
---|
1397 | + ( w_comp - ( w(k-1,j-1,i) + w(k-1,j,i) ) ) * ddzw(k) & |
---|
1398 | ) * 0.5 |
---|
1399 | |
---|
1400 | tend(k,j,i) = tend(k,j,i) - ( & |
---|
1401 | ( flux_r(k) + diss_r(k) & |
---|
1402 | - flux_l_v(k,j,tn) - diss_l_v(k,j,tn) ) * ddx & |
---|
1403 | + ( flux_n(k) + diss_n(k) & |
---|
1404 | - flux_s_v(k,tn) - diss_s_v(k,tn) ) * ddy & |
---|
1405 | + ( flux_t(k) + diss_t(k) & |
---|
1406 | - flux_d - diss_d ) * ddzw(k) & |
---|
1407 | ) + v(k,j,i) * div |
---|
1408 | |
---|
1409 | flux_l_v(k,j,tn) = flux_r(k) |
---|
1410 | diss_l_v(k,j,tn) = diss_r(k) |
---|
1411 | flux_s_v(k,tn) = flux_n(k) |
---|
1412 | diss_s_v(k,tn) = diss_n(k) |
---|
1413 | flux_d = flux_t(k) |
---|
1414 | diss_d = diss_t(k) |
---|
1415 | |
---|
1416 | ! |
---|
1417 | !-- Statistical Evaluation of v'v'. The factor has to be applied for |
---|
1418 | !-- right evaluation when gallilei_trans = .T. . |
---|
1419 | sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) & |
---|
1420 | + ( flux_n(k) & |
---|
1421 | * ( v_comp(k) - 2.0 * hom(k,1,2,0) ) & |
---|
1422 | / ( v_comp(k) - gv + 1.0E-20 ) & |
---|
1423 | + diss_n(k) & |
---|
1424 | * ABS( v_comp(k) - 2.0 * hom(k,1,2,0) ) & |
---|
1425 | / ( ABS( v_comp(k) - gv ) +1.0E-20 ) ) & |
---|
1426 | * weight_substep(intermediate_timestep_count) |
---|
1427 | ! |
---|
1428 | !-- Statistical Evaluation of w'v'. |
---|
1429 | sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) & |
---|
1430 | + ( flux_t(k) + diss_t(k) ) & |
---|
1431 | * weight_substep(intermediate_timestep_count) |
---|
1432 | |
---|
1433 | ENDDO |
---|
1434 | sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn) |
---|
1435 | |
---|
1436 | |
---|
1437 | END SUBROUTINE advec_v_ws_ij |
---|
1438 | |
---|
1439 | |
---|
1440 | |
---|
1441 | !------------------------------------------------------------------------------! |
---|
1442 | ! Advection of w-component - Call for grid point i,j |
---|
1443 | !------------------------------------------------------------------------------! |
---|
1444 | SUBROUTINE advec_w_ws_ij( i, j, i_omp, tn ) |
---|
1445 | |
---|
1446 | USE arrays_3d |
---|
1447 | USE constants |
---|
1448 | USE control_parameters |
---|
1449 | USE grid_variables |
---|
1450 | USE indices |
---|
1451 | USE statistics |
---|
1452 | |
---|
1453 | IMPLICIT NONE |
---|
1454 | |
---|
1455 | INTEGER :: i, i_omp, j, k, k_mm, k_pp, k_ppp, tn |
---|
1456 | REAL :: diss_d, div, flux_d, gu, gv, u_comp, v_comp, w_comp |
---|
1457 | REAL, DIMENSION(nzb:nzt+1) :: diss_n, diss_r, diss_t, flux_n, flux_r, & |
---|
1458 | flux_t |
---|
1459 | |
---|
1460 | gu = 2.0 * u_gtrans |
---|
1461 | gv = 2.0 * v_gtrans |
---|
1462 | |
---|
1463 | ! |
---|
1464 | !-- Compute southside fluxes for the respective boundary. |
---|
1465 | IF ( j == nys ) THEN |
---|
1466 | |
---|
1467 | DO k = nzb+1, nzb_max |
---|
1468 | |
---|
1469 | v_comp = v(k+1,j,i) + v(k,j,i) - gv |
---|
1470 | flux_s_w(k,tn) = v_comp * ( & |
---|
1471 | ( 37.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 & |
---|
1472 | + 7.0 * IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3 & |
---|
1473 | + IBITS(wall_flags_0(k,j,i),30,1) * adv_mom_1 & |
---|
1474 | ) * & |
---|
1475 | ( w(k,j,i) + w(k,j-1,i) ) & |
---|
1476 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 & |
---|
1477 | + IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3 & |
---|
1478 | ) * & |
---|
1479 | ( w(k,j+1,i) + w(k,j-2,i) ) & |
---|
1480 | + ( IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 & |
---|
1481 | ) * & |
---|
1482 | ( w(k,j+2,i) + w(k,j-3,i) ) & |
---|
1483 | ) |
---|
1484 | |
---|
1485 | diss_s_w(k,tn) = - ABS( v_comp ) * ( & |
---|
1486 | ( 10.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 & |
---|
1487 | + 3.0 * IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3 & |
---|
1488 | + IBITS(wall_flags_0(k,j,i),30,1) * adv_mom_1 & |
---|
1489 | ) * & |
---|
1490 | ( w(k,j,i) - w(k,j-1,i) ) & |
---|
1491 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 & |
---|
1492 | + IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3 & |
---|
1493 | ) * & |
---|
1494 | ( w(k,j+1,i) - w(k,j-2,i) ) & |
---|
1495 | + ( IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 & |
---|
1496 | ) * & |
---|
1497 | ( w(k,j+2,i) - w(k,j-3,i) ) & |
---|
1498 | ) |
---|
1499 | |
---|
1500 | ENDDO |
---|
1501 | |
---|
1502 | DO k = nzb_max+1, nzt |
---|
1503 | |
---|
1504 | v_comp = v(k+1,j,i) + v(k,j,i) - gv |
---|
1505 | flux_s_w(k,tn) = v_comp * ( & |
---|
1506 | 37.0 * ( w(k,j,i) + w(k,j-1,i) ) & |
---|
1507 | - 8.0 * ( w(k,j+1,i) +w(k,j-2,i) ) & |
---|
1508 | + ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5 |
---|
1509 | diss_s_w(k,tn) = - ABS( v_comp ) * ( & |
---|
1510 | 10.0 * ( w(k,j,i) - w(k,j-1,i) ) & |
---|
1511 | - 5.0 * ( w(k,j+1,i) - w(k,j-2,i) ) & |
---|
1512 | + ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5 |
---|
1513 | |
---|
1514 | ENDDO |
---|
1515 | |
---|
1516 | ENDIF |
---|
1517 | ! |
---|
1518 | !-- Compute leftside fluxes for the respective boundary. |
---|
1519 | IF ( i == i_omp ) THEN |
---|
1520 | |
---|
1521 | DO k = nzb+1, nzb_max |
---|
1522 | |
---|
1523 | u_comp = u(k+1,j,i) + u(k,j,i) - gu |
---|
1524 | flux_l_w(k,j,tn) = u_comp * ( & |
---|
1525 | ( 37.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 & |
---|
1526 | + 7.0 * IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3 & |
---|
1527 | + IBITS(wall_flags_0(k,j,i),27,1) * adv_mom_1 & |
---|
1528 | ) * & |
---|
1529 | ( w(k,j,i) + w(k,j,i-1) ) & |
---|
1530 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 & |
---|
1531 | + IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3 & |
---|
1532 | ) * & |
---|
1533 | ( w(k,j,i+1) + w(k,j,i-2) ) & |
---|
1534 | + ( IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 & |
---|
1535 | ) * & |
---|
1536 | ( w(k,j,i+2) + w(k,j,i-3) ) & |
---|
1537 | ) |
---|
1538 | |
---|
1539 | diss_l_w(k,j,tn) = - ABS( u_comp ) * ( & |
---|
1540 | ( 10.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 & |
---|
1541 | + 3.0 * IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3 & |
---|
1542 | + IBITS(wall_flags_0(k,j,i),27,1) * adv_mom_1 & |
---|
1543 | ) * & |
---|
1544 | ( w(k,j,i) - w(k,j,i-1) ) & |
---|
1545 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 & |
---|
1546 | + IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3 & |
---|
1547 | ) * & |
---|
1548 | ( w(k,j,i+1) - w(k,j,i-2) ) & |
---|
1549 | + ( IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 & |
---|
1550 | ) * & |
---|
1551 | ( w(k,j,i+2) - w(k,j,i-3) ) & |
---|
1552 | ) |
---|
1553 | |
---|
1554 | ENDDO |
---|
1555 | |
---|
1556 | DO k = nzb_max+1, nzt |
---|
1557 | |
---|
1558 | u_comp = u(k+1,j,i) + u(k,j,i) - gu |
---|
1559 | flux_l_w(k,j,tn) = u_comp * ( & |
---|
1560 | 37.0 * ( w(k,j,i) + w(k,j,i-1) ) & |
---|
1561 | - 8.0 * ( w(k,j,i+1) + w(k,j,i-2) ) & |
---|
1562 | + ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5 |
---|
1563 | diss_l_w(k,j,tn) = - ABS( u_comp ) * ( & |
---|
1564 | 10.0 * ( w(k,j,i) - w(k,j,i-1) ) & |
---|
1565 | - 5.0 * ( w(k,j,i+1) - w(k,j,i-2) ) & |
---|
1566 | + ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5 |
---|
1567 | |
---|
1568 | ENDDO |
---|
1569 | |
---|
1570 | ENDIF |
---|
1571 | ! |
---|
1572 | !-- The lower flux has to be calculated explicetely for the tendency at |
---|
1573 | !-- the first w-level. For topography wall this is done implicitely by |
---|
1574 | !-- wall_flags_0. |
---|
1575 | k = nzb + 1 |
---|
1576 | w_comp = w(k,j,i) + w(k-1,j,i) |
---|
1577 | flux_t(0) = w_comp * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1 |
---|
1578 | diss_t(0) = -ABS(w_comp) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1 |
---|
1579 | flux_d = flux_t(0) |
---|
1580 | diss_d = diss_t(0) |
---|
1581 | ! |
---|
1582 | !-- Now compute the fluxes and tendency terms for the horizontal |
---|
1583 | !-- and vertical parts. |
---|
1584 | DO k = nzb+1, nzb_max |
---|
1585 | |
---|
1586 | u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu |
---|
1587 | flux_r(k) = u_comp * ( & |
---|
1588 | ( 37.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 & |
---|
1589 | + 7.0 * IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3 & |
---|
1590 | + IBITS(wall_flags_0(k,j,i),27,1) * adv_mom_1 & |
---|
1591 | ) * & |
---|
1592 | ( w(k,j,i+1) + w(k,j,i) ) & |
---|
1593 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 & |
---|
1594 | + IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3 & |
---|
1595 | ) * & |
---|
1596 | ( w(k,j,i+2) + w(k,j,i-1) ) & |
---|
1597 | + ( IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 & |
---|
1598 | ) * & |
---|
1599 | ( w(k,j,i+3) + w(k,j,i-2) ) & |
---|
1600 | ) |
---|
1601 | |
---|
1602 | diss_r(k) = - ABS( u_comp ) * ( & |
---|
1603 | ( 10.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 & |
---|
1604 | + 3.0 * IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3 & |
---|
1605 | + IBITS(wall_flags_0(k,j,i),27,1) * adv_mom_1 & |
---|
1606 | ) * & |
---|
1607 | ( w(k,j,i+1) - w(k,j,i) ) & |
---|
1608 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 & |
---|
1609 | + IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3 & |
---|
1610 | ) * & |
---|
1611 | ( w(k,j,i+2) - w(k,j,i-1) ) & |
---|
1612 | + ( IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 & |
---|
1613 | ) * & |
---|
1614 | ( w(k,j,i+3) - w(k,j,i-2) ) & |
---|
1615 | ) |
---|
1616 | |
---|
1617 | v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv |
---|
1618 | flux_n(k) = v_comp * ( & |
---|
1619 | ( 37.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 & |
---|
1620 | + 7.0 * IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3 & |
---|
1621 | + IBITS(wall_flags_0(k,j,i),30,1) * adv_mom_1 & |
---|
1622 | ) * & |
---|
1623 | ( w(k,j+1,i) + w(k,j,i) ) & |
---|
1624 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 & |
---|
1625 | + IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3 & |
---|
1626 | ) * & |
---|
1627 | ( w(k,j+2,i) + w(k,j-1,i) ) & |
---|
1628 | + ( IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 & |
---|
1629 | ) * & |
---|
1630 | ( w(k,j+3,i) + w(k,j-2,i) ) & |
---|
1631 | ) |
---|
1632 | |
---|
1633 | diss_n(k) = - ABS( v_comp ) * ( & |
---|
1634 | ( 10.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 & |
---|
1635 | + 3.0 * IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3 & |
---|
1636 | + IBITS(wall_flags_0(k,j,i),30,1) * adv_mom_1 & |
---|
1637 | ) * & |
---|
1638 | ( w(k,j+1,i) - w(k,j,i) ) & |
---|
1639 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 & |
---|
1640 | + IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3 & |
---|
1641 | ) * & |
---|
1642 | ( w(k,j+2,i) - w(k,j-1,i) ) & |
---|
1643 | + ( IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 & |
---|
1644 | ) * & |
---|
1645 | ( w(k,j+3,i) - w(k,j-2,i) ) & |
---|
1646 | ) |
---|
1647 | ! |
---|
1648 | !-- k index has to be modified near bottom and top, else array |
---|
1649 | !-- subscripts will be exceeded. |
---|
1650 | k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),35,1) |
---|
1651 | k_pp = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),33,1) ) |
---|
1652 | k_mm = k - 2 * IBITS(wall_flags_0(k,j,i),35,1) |
---|
1653 | |
---|
1654 | w_comp = w(k+1,j,i) + w(k,j,i) |
---|
1655 | flux_t(k) = w_comp * ( & |
---|
1656 | ( 37.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5 & |
---|
1657 | + 7.0 * IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3 & |
---|
1658 | + IBITS(wall_flags_0(k,j,i),33,1) * adv_mom_1 & |
---|
1659 | ) * & |
---|
1660 | ( w(k+1,j,i) + w(k,j,i) ) & |
---|
1661 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5 & |
---|
1662 | + IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3 & |
---|
1663 | ) * & |
---|
1664 | ( w(k_pp,j,i) + w(k-1,j,i) ) & |
---|
1665 | + ( IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5 & |
---|
1666 | ) * & |
---|
1667 | ( w(k_ppp,j,i) + w(k_mm,j,i) ) & |
---|
1668 | ) |
---|
1669 | |
---|
1670 | diss_t(k) = - ABS( w_comp ) * ( & |
---|
1671 | ( 10.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5 & |
---|
1672 | + 3.0 * IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3 & |
---|
1673 | + IBITS(wall_flags_0(k,j,i),33,1) * adv_mom_1 & |
---|
1674 | ) * & |
---|
1675 | ( w(k+1,j,i) - w(k,j,i) ) & |
---|
1676 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5 & |
---|
1677 | + IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3 & |
---|
1678 | ) * & |
---|
1679 | ( w(k_pp,j,i) - w(k-1,j,i) ) & |
---|
1680 | + ( IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5 & |
---|
1681 | ) * & |
---|
1682 | ( w(k_ppp,j,i) - w(k_mm,j,i) ) & |
---|
1683 | ) |
---|
1684 | |
---|
1685 | ! |
---|
1686 | !-- Calculate the divergence of the velocity field. A respective |
---|
1687 | !-- correction is needed to overcome numerical instabilities introduced |
---|
1688 | !-- by a not sufficient reduction of divergences near topography. |
---|
1689 | div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx & |
---|
1690 | + ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy & |
---|
1691 | + ( w_comp - ( w(k,j,i) + w(k-1,j,i) ) ) * ddzu(k+1) & |
---|
1692 | ) * 0.5 |
---|
1693 | |
---|
1694 | tend(k,j,i) = tend(k,j,i) - ( & |
---|
1695 | ( flux_r(k) + diss_r(k) & |
---|
1696 | - flux_l_w(k,j,tn) - diss_l_w(k,j,tn) ) * ddx & |
---|
1697 | + ( flux_n(k) + diss_n(k) & |
---|
1698 | - flux_s_w(k,tn) - diss_s_w(k,tn) ) * ddy & |
---|
1699 | + ( flux_t(k) + diss_t(k) & |
---|
1700 | - flux_d - diss_d ) * ddzu(k+1) & |
---|
1701 | ) + div * w(k,j,i) |
---|
1702 | |
---|
1703 | flux_l_w(k,j,tn) = flux_r(k) |
---|
1704 | diss_l_w(k,j,tn) = diss_r(k) |
---|
1705 | flux_s_w(k,tn) = flux_n(k) |
---|
1706 | diss_s_w(k,tn) = diss_n(k) |
---|
1707 | flux_d = flux_t(k) |
---|
1708 | diss_d = diss_t(k) |
---|
1709 | ! |
---|
1710 | !-- Statistical Evaluation of w'w'. |
---|
1711 | sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) & |
---|
1712 | + ( flux_t(k) + diss_t(k) ) & |
---|
1713 | * weight_substep(intermediate_timestep_count) |
---|
1714 | |
---|
1715 | ENDDO |
---|
1716 | |
---|
1717 | DO k = nzb_max+1, nzt |
---|
1718 | |
---|
1719 | u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu |
---|
1720 | flux_r(k) = u_comp * ( & |
---|
1721 | 37.0 * ( w(k,j,i+1) + w(k,j,i) ) & |
---|
1722 | - 8.0 * ( w(k,j,i+2) + w(k,j,i-1) ) & |
---|
1723 | + ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5 |
---|
1724 | |
---|
1725 | diss_r(k) = - ABS( u_comp ) * ( & |
---|
1726 | 10.0 * ( w(k,j,i+1) - w(k,j,i) ) & |
---|
1727 | - 5.0 * ( w(k,j,i+2) - w(k,j,i-1) ) & |
---|
1728 | + ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5 |
---|
1729 | |
---|
1730 | v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv |
---|
1731 | flux_n(k) = v_comp * ( & |
---|
1732 | 37.0 * ( w(k,j+1,i) + w(k,j,i) ) & |
---|
1733 | - 8.0 * ( w(k,j+2,i) + w(k,j-1,i) ) & |
---|
1734 | + ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5 |
---|
1735 | |
---|
1736 | diss_n(k) = - ABS( v_comp ) * ( & |
---|
1737 | 10.0 * ( w(k,j+1,i) - w(k,j,i) ) & |
---|
1738 | - 5.0 * ( w(k,j+2,i) - w(k,j-1,i) ) & |
---|
1739 | + ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5 |
---|
1740 | ! |
---|
1741 | !-- k index has to be modified near bottom and top, else array |
---|
1742 | !-- subscripts will be exceeded. |
---|
1743 | k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),35,1) |
---|
1744 | k_pp = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),33,1) ) |
---|
1745 | k_mm = k - 2 * IBITS(wall_flags_0(k,j,i),35,1) |
---|
1746 | |
---|
1747 | w_comp = w(k+1,j,i) + w(k,j,i) |
---|
1748 | flux_t(k) = w_comp * ( & |
---|
1749 | ( 37.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5 & |
---|
1750 | + 7.0 * IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3 & |
---|
1751 | + IBITS(wall_flags_0(k,j,i),33,1) * adv_mom_1 & |
---|
1752 | ) * & |
---|
1753 | ( w(k+1,j,i) + w(k,j,i) ) & |
---|
1754 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5 & |
---|
1755 | + IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3 & |
---|
1756 | ) * & |
---|
1757 | ( w(k_pp,j,i) + w(k-1,j,i) ) & |
---|
1758 | + ( IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5 & |
---|
1759 | ) * & |
---|
1760 | ( w(k_ppp,j,i) + w(k_mm,j,i) ) & |
---|
1761 | ) |
---|
1762 | |
---|
1763 | diss_t(k) = - ABS( w_comp ) * ( & |
---|
1764 | ( 10.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5 & |
---|
1765 | + 3.0 * IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3 & |
---|
1766 | + IBITS(wall_flags_0(k,j,i),33,1) * adv_mom_1 & |
---|
1767 | ) * & |
---|
1768 | ( w(k+1,j,i) - w(k,j,i) ) & |
---|
1769 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5 & |
---|
1770 | + IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3 & |
---|
1771 | ) * & |
---|
1772 | ( w(k_pp,j,i) - w(k-1,j,i) ) & |
---|
1773 | + ( IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5 & |
---|
1774 | ) * & |
---|
1775 | ( w(k_ppp,j,i) - w(k_mm,j,i) ) & |
---|
1776 | ) |
---|
1777 | ! |
---|
1778 | !-- Calculate the divergence of the velocity field. A respective |
---|
1779 | !-- correction is needed to overcome numerical instabilities introduced |
---|
1780 | !-- by a not sufficient reduction of divergences near topography. |
---|
1781 | div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx & |
---|
1782 | + ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy & |
---|
1783 | + ( w_comp - ( w(k,j,i) + w(k-1,j,i) ) ) * ddzu(k+1) & |
---|
1784 | ) * 0.5 |
---|
1785 | |
---|
1786 | tend(k,j,i) = tend(k,j,i) - ( & |
---|
1787 | ( flux_r(k) + diss_r(k) & |
---|
1788 | - flux_l_w(k,j,tn) - diss_l_w(k,j,tn) ) * ddx & |
---|
1789 | + ( flux_n(k) + diss_n(k) & |
---|
1790 | - flux_s_w(k,tn) - diss_s_w(k,tn) ) * ddy & |
---|
1791 | + ( flux_t(k) + diss_t(k) & |
---|
1792 | - flux_d - diss_d ) * ddzu(k+1) & |
---|
1793 | ) + div * w(k,j,i) |
---|
1794 | |
---|
1795 | flux_l_w(k,j,tn) = flux_r(k) |
---|
1796 | diss_l_w(k,j,tn) = diss_r(k) |
---|
1797 | flux_s_w(k,tn) = flux_n(k) |
---|
1798 | diss_s_w(k,tn) = diss_n(k) |
---|
1799 | flux_d = flux_t(k) |
---|
1800 | diss_d = diss_t(k) |
---|
1801 | ! |
---|
1802 | !-- Statistical Evaluation of w'w'. |
---|
1803 | sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) & |
---|
1804 | + ( flux_t(k) + diss_t(k) ) & |
---|
1805 | * weight_substep(intermediate_timestep_count) |
---|
1806 | |
---|
1807 | ENDDO |
---|
1808 | |
---|
1809 | |
---|
1810 | END SUBROUTINE advec_w_ws_ij |
---|
1811 | |
---|
1812 | |
---|
1813 | !------------------------------------------------------------------------------! |
---|
1814 | ! Scalar advection - Call for all grid points |
---|
1815 | !------------------------------------------------------------------------------! |
---|
1816 | SUBROUTINE advec_s_ws( sk, sk_char ) |
---|
1817 | |
---|
1818 | USE arrays_3d |
---|
1819 | USE constants |
---|
1820 | USE control_parameters |
---|
1821 | USE grid_variables |
---|
1822 | USE indices |
---|
1823 | USE statistics |
---|
1824 | |
---|
1825 | IMPLICIT NONE |
---|
1826 | |
---|
1827 | INTEGER :: i, j, k, k_mm, k_pp, k_ppp, tn = 0 |
---|
1828 | REAL, DIMENSION(:,:,:), POINTER :: sk |
---|
1829 | REAL :: diss_d, div, flux_d, u_comp, v_comp |
---|
1830 | REAL, DIMENSION(nzb:nzt) :: diss_n, diss_r, diss_t, flux_n, flux_r, & |
---|
1831 | flux_t |
---|
1832 | REAL, DIMENSION(nzb+1:nzt) :: swap_diss_y_local, swap_flux_y_local |
---|
1833 | REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local, & |
---|
1834 | swap_flux_x_local |
---|
1835 | CHARACTER (LEN = *), INTENT(IN) :: sk_char |
---|
1836 | |
---|
1837 | ! |
---|
1838 | !-- Compute the fluxes for the whole left boundary of the processor domain. |
---|
1839 | i = nxl |
---|
1840 | DO j = nys, nyn |
---|
1841 | |
---|
1842 | DO k = nzb+1, nzb_max |
---|
1843 | |
---|
1844 | u_comp = u(k,j,i) - u_gtrans |
---|
1845 | swap_flux_x_local(k,j) = u_comp * ( & |
---|
1846 | ( 37.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 & |
---|
1847 | + 7.0 * IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3 & |
---|
1848 | + IBITS(wall_flags_0(k,j,i),0,1) * adv_sca_1 & |
---|
1849 | ) * & |
---|
1850 | ( sk(k,j,i) + sk(k,j,i-1) ) & |
---|
1851 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 & |
---|
1852 | + IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3 & |
---|
1853 | ) * & |
---|
1854 | ( sk(k,j,i+1) + sk(k,j,i-2) ) & |
---|
1855 | + ( IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 & |
---|
1856 | ) * ( sk(k,j,i+2) + sk(k,j,i-3) ) & |
---|
1857 | ) |
---|
1858 | |
---|
1859 | swap_diss_x_local(k,j) = -ABS( u_comp ) * ( & |
---|
1860 | ( 10.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 & |
---|
1861 | + 3.0 * IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3 & |
---|
1862 | + IBITS(wall_flags_0(k,j,i),0,1) * adv_sca_1 & |
---|
1863 | ) * & |
---|
1864 | ( sk(k,j,i) - sk(k,j,i-1) ) & |
---|
1865 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 & |
---|
1866 | + IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3 & |
---|
1867 | ) * & |
---|
1868 | ( sk(k,j,i+1) - sk(k,j,i-2) ) & |
---|
1869 | + ( IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 & |
---|
1870 | ) * & |
---|
1871 | ( sk(k,j,i+2) - sk(k,j,i-3) ) & |
---|
1872 | ) |
---|
1873 | |
---|
1874 | ENDDO |
---|
1875 | |
---|
1876 | DO k = nzb_max+1, nzt |
---|
1877 | |
---|
1878 | u_comp = u(k,j,i) - u_gtrans |
---|
1879 | swap_flux_x_local(k,j) = u_comp * ( & |
---|
1880 | 37.0 * ( sk(k,j,i) + sk(k,j,i-1) ) & |
---|
1881 | - 8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) ) & |
---|
1882 | + ( sk(k,j,i+2) + sk(k,j,i-3) ) & |
---|
1883 | ) * adv_sca_5 |
---|
1884 | |
---|
1885 | swap_diss_x_local(k,j) = -ABS( u_comp ) * ( & |
---|
1886 | 10.0 * ( sk(k,j,i) - sk(k,j,i-1) ) & |
---|
1887 | - 5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) ) & |
---|
1888 | + ( sk(k,j,i+2) - sk(k,j,i-3) ) & |
---|
1889 | ) * adv_sca_5 |
---|
1890 | |
---|
1891 | ENDDO |
---|
1892 | |
---|
1893 | ENDDO |
---|
1894 | |
---|
1895 | DO i = nxl, nxr |
---|
1896 | |
---|
1897 | j = nys |
---|
1898 | DO k = nzb+1, nzb_max |
---|
1899 | |
---|
1900 | v_comp = v(k,j,i) - v_gtrans |
---|
1901 | swap_flux_y_local(k) = v_comp * ( & |
---|
1902 | ( 37.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 & |
---|
1903 | + 7.0 * IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3 & |
---|
1904 | + IBITS(wall_flags_0(k,j,i),3,1) * adv_sca_1 & |
---|
1905 | ) * & |
---|
1906 | ( sk(k,j,i) + sk(k,j-1,i) ) & |
---|
1907 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 & |
---|
1908 | + IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3 & |
---|
1909 | ) * & |
---|
1910 | ( sk(k,j+1,i) + sk(k,j-2,i) ) & |
---|
1911 | + ( IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 & |
---|
1912 | ) * & |
---|
1913 | ( sk(k,j+2,i) + sk(k,j-3,i) ) & |
---|
1914 | ) |
---|
1915 | |
---|
1916 | swap_diss_y_local(k) = -ABS( v_comp ) * ( & |
---|
1917 | ( 10.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 & |
---|
1918 | + 3.0 * IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3 & |
---|
1919 | + IBITS(wall_flags_0(k,j,i),3,1) * adv_sca_1 & |
---|
1920 | ) * & |
---|
1921 | ( sk(k,j,i) - sk(k,j-1,i) ) & |
---|
1922 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 & |
---|
1923 | + IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3 & |
---|
1924 | ) * & |
---|
1925 | ( sk(k,j+1,i) - sk(k,j-2,i) ) & |
---|
1926 | + ( IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 & |
---|
1927 | ) * & |
---|
1928 | ( sk(k,j+2,i) - sk(k,j-3,i) ) & |
---|
1929 | ) |
---|
1930 | |
---|
1931 | ENDDO |
---|
1932 | ! |
---|
1933 | !-- Above to the top of the highest topography. No degradation necessary. |
---|
1934 | DO k = nzb_max+1, nzt |
---|
1935 | |
---|
1936 | v_comp = v(k,j,i) - v_gtrans |
---|
1937 | swap_flux_y_local(k) = v_comp * ( & |
---|
1938 | 37.0 * ( sk(k,j,i) + sk(k,j-1,i) ) & |
---|
1939 | - 8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) ) & |
---|
1940 | + ( sk(k,j+2,i) + sk(k,j-3,i) ) & |
---|
1941 | ) * adv_sca_5 |
---|
1942 | swap_diss_y_local(k) = -ABS( v_comp ) * ( & |
---|
1943 | 10.0 * ( sk(k,j,i) - sk(k,j-1,i) ) & |
---|
1944 | - 5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) ) & |
---|
1945 | + sk(k,j+2,i) - sk(k,j-3,i) & |
---|
1946 | ) * adv_sca_5 |
---|
1947 | |
---|
1948 | ENDDO |
---|
1949 | |
---|
1950 | DO j = nys, nyn |
---|
1951 | |
---|
1952 | flux_t(0) = 0.0 |
---|
1953 | diss_t(0) = 0.0 |
---|
1954 | flux_d = 0.0 |
---|
1955 | diss_d = 0.0 |
---|
1956 | |
---|
1957 | DO k = nzb+1, nzb_max |
---|
1958 | |
---|
1959 | u_comp = u(k,j,i+1) - u_gtrans |
---|
1960 | flux_r(k) = u_comp * ( & |
---|
1961 | ( 37.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 & |
---|
1962 | + 7.0 * IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3 & |
---|
1963 | + IBITS(wall_flags_0(k,j,i),0,1) * adv_sca_1 & |
---|
1964 | ) * & |
---|
1965 | ( sk(k,j,i+1) + sk(k,j,i) ) & |
---|
1966 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 & |
---|
1967 | + IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3 & |
---|
1968 | ) * & |
---|
1969 | ( sk(k,j,i+2) + sk(k,j,i-1) ) & |
---|
1970 | + ( IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 & |
---|
1971 | ) * ( sk(k,j,i+3) + sk(k,j,i-2) ) & |
---|
1972 | ) |
---|
1973 | |
---|
1974 | diss_r(k) = -ABS( u_comp ) * ( & |
---|
1975 | ( 10.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 & |
---|
1976 | + 3.0 * IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3 & |
---|
1977 | + IBITS(wall_flags_0(k,j,i),0,1) * adv_sca_1 & |
---|
1978 | ) * & |
---|
1979 | ( sk(k,j,i+1) - sk(k,j,i) ) & |
---|
1980 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 & |
---|
1981 | + IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3 & |
---|
1982 | ) * & |
---|
1983 | ( sk(k,j,i+2) - sk(k,j,i-1) ) & |
---|
1984 | + ( IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 & |
---|
1985 | ) * & |
---|
1986 | ( sk(k,j,i+3) - sk(k,j,i-2) ) & |
---|
1987 | ) |
---|
1988 | |
---|
1989 | v_comp = v(k,j+1,i) - v_gtrans |
---|
1990 | flux_n(k) = v_comp * ( & |
---|
1991 | ( 37.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 & |
---|
1992 | + 7.0 * IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3 & |
---|
1993 | + IBITS(wall_flags_0(k,j,i),3,1) * adv_sca_1 & |
---|
1994 | ) * & |
---|
1995 | ( sk(k,j+1,i) + sk(k,j,i) ) & |
---|
1996 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 & |
---|
1997 | + IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3 & |
---|
1998 | ) * & |
---|
1999 | ( sk(k,j+2,i) + sk(k,j-1,i) ) & |
---|
2000 | + ( IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 & |
---|
2001 | ) * & |
---|
2002 | ( sk(k,j+3,i) + sk(k,j-2,i) ) & |
---|
2003 | ) |
---|
2004 | |
---|
2005 | diss_n(k) = -ABS( v_comp ) * ( & |
---|
2006 | ( 10.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 & |
---|
2007 | + 3.0 * IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3 & |
---|
2008 | + IBITS(wall_flags_0(k,j,i),3,1) * adv_sca_1 & |
---|
2009 | ) * & |
---|
2010 | ( sk(k,j+1,i) - sk(k,j,i) ) & |
---|
2011 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 & |
---|
2012 | + IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3 & |
---|
2013 | ) * & |
---|
2014 | ( sk(k,j+2,i) - sk(k,j-1,i) ) & |
---|
2015 | + ( IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 & |
---|
2016 | ) * & |
---|
2017 | ( sk(k,j+3,i) - sk(k,j-2,i) ) & |
---|
2018 | ) |
---|
2019 | ! |
---|
2020 | !-- k index has to be modified near bottom and top, else array |
---|
2021 | !-- subscripts will be exceeded. |
---|
2022 | k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),8,1) |
---|
2023 | k_pp = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),6,1) ) |
---|
2024 | k_mm = k - 2 * IBITS(wall_flags_0(k,j,i),8,1) |
---|
2025 | |
---|
2026 | |
---|
2027 | flux_t(k) = w(k,j,i) * ( & |
---|
2028 | ( 37.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 & |
---|
2029 | + 7.0 * IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3 & |
---|
2030 | + IBITS(wall_flags_0(k,j,i),6,1) * adv_sca_1 & |
---|
2031 | ) * & |
---|
2032 | ( sk(k+1,j,i) + sk(k,j,i) ) & |
---|
2033 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 & |
---|
2034 | + IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3 & |
---|
2035 | ) * & |
---|
2036 | ( sk(k_pp,j,i) + sk(k-1,j,i) ) & |
---|
2037 | + ( IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 & |
---|
2038 | ) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) & |
---|
2039 | ) |
---|
2040 | |
---|
2041 | diss_t(k) = -ABS( w(k,j,i) ) * ( & |
---|
2042 | ( 10.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 & |
---|
2043 | + 3.0 * IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3 & |
---|
2044 | + IBITS(wall_flags_0(k,j,i),6,1) * adv_sca_1 & |
---|
2045 | ) * & |
---|
2046 | ( sk(k+1,j,i) - sk(k,j,i) ) & |
---|
2047 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 & |
---|
2048 | + IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3 & |
---|
2049 | ) * & |
---|
2050 | ( sk(k_pp,j,i) - sk(k-1,j,i) ) & |
---|
2051 | + ( IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 & |
---|
2052 | ) * & |
---|
2053 | ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & |
---|
2054 | ) |
---|
2055 | ! |
---|
2056 | !-- Calculate the divergence of the velocity field. A respective |
---|
2057 | !-- correction is needed to overcome numerical instabilities caused |
---|
2058 | !-- by a not sufficient reduction of divergences near topography. |
---|
2059 | div = ( u(k,j,i+1) - u(k,j,i) ) * ddx & |
---|
2060 | + ( v(k,j+1,i) - v(k,j,i) ) * ddy & |
---|
2061 | + ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) |
---|
2062 | |
---|
2063 | tend(k,j,i) = tend(k,j,i) - ( & |
---|
2064 | ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) - & |
---|
2065 | swap_diss_x_local(k,j) ) * ddx & |
---|
2066 | + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k) - & |
---|
2067 | swap_diss_y_local(k) ) * ddy & |
---|
2068 | + ( flux_t(k) + diss_t(k) - flux_d - diss_d & |
---|
2069 | ) * ddzw(k) & |
---|
2070 | ) + sk(k,j,i) * div |
---|
2071 | |
---|
2072 | swap_flux_y_local(k) = flux_n(k) |
---|
2073 | swap_diss_y_local(k) = diss_n(k) |
---|
2074 | swap_flux_x_local(k,j) = flux_r(k) |
---|
2075 | swap_diss_x_local(k,j) = diss_r(k) |
---|
2076 | flux_d = flux_t(k) |
---|
2077 | diss_d = diss_t(k) |
---|
2078 | |
---|
2079 | ENDDO |
---|
2080 | |
---|
2081 | DO k = nzb_max+1, nzt |
---|
2082 | |
---|
2083 | u_comp = u(k,j,i+1) - u_gtrans |
---|
2084 | flux_r(k) = u_comp * ( & |
---|
2085 | 37.0 * ( sk(k,j,i+1) + sk(k,j,i) ) & |
---|
2086 | - 8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) ) & |
---|
2087 | + ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5 |
---|
2088 | diss_r(k) = -ABS( u_comp ) * ( & |
---|
2089 | 10.0 * ( sk(k,j,i+1) - sk(k,j,i) ) & |
---|
2090 | - 5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) ) & |
---|
2091 | + ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5 |
---|
2092 | |
---|
2093 | v_comp = v(k,j+1,i) - v_gtrans |
---|
2094 | flux_n(k) = v_comp * ( & |
---|
2095 | 37.0 * ( sk(k,j+1,i) + sk(k,j,i) ) & |
---|
2096 | - 8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) ) & |
---|
2097 | + ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5 |
---|
2098 | diss_n(k) = -ABS( v_comp ) * ( & |
---|
2099 | 10.0 * ( sk(k,j+1,i) - sk(k,j,i) ) & |
---|
2100 | - 5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) ) & |
---|
2101 | + ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5 |
---|
2102 | ! |
---|
2103 | !-- k index has to be modified near bottom and top, else array |
---|
2104 | !-- subscripts will be exceeded. |
---|
2105 | k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),8,1) |
---|
2106 | k_pp = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),6,1) ) |
---|
2107 | k_mm = k - 2 * IBITS(wall_flags_0(k,j,i),8,1) |
---|
2108 | |
---|
2109 | flux_t(k) = w(k,j,i) * ( & |
---|
2110 | ( 37.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 & |
---|
2111 | + 7.0 * IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3 & |
---|
2112 | + IBITS(wall_flags_0(k,j,i),6,1) * adv_sca_1 & |
---|
2113 | ) * & |
---|
2114 | ( sk(k+1,j,i) + sk(k,j,i) ) & |
---|
2115 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 & |
---|
2116 | + IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3 & |
---|
2117 | ) * & |
---|
2118 | ( sk(k_pp,j,i) + sk(k-1,j,i) ) & |
---|
2119 | + ( IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 & |
---|
2120 | ) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) & |
---|
2121 | ) |
---|
2122 | |
---|
2123 | diss_t(k) = -ABS( w(k,j,i) ) * ( & |
---|
2124 | ( 10.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 & |
---|
2125 | + 3.0 * IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3 & |
---|
2126 | + IBITS(wall_flags_0(k,j,i),6,1) * adv_sca_1 & |
---|
2127 | ) * & |
---|
2128 | ( sk(k+1,j,i) - sk(k,j,i) ) & |
---|
2129 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 & |
---|
2130 | + IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3 & |
---|
2131 | ) * & |
---|
2132 | ( sk(k_pp,j,i) - sk(k-1,j,i) ) & |
---|
2133 | + ( IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 & |
---|
2134 | ) * & |
---|
2135 | ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & |
---|
2136 | ) |
---|
2137 | ! |
---|
2138 | !-- Calculate the divergence of the velocity field. A respective |
---|
2139 | !-- correction is needed to overcome numerical instabilities introduced |
---|
2140 | !-- by a not sufficient reduction of divergences near topography. |
---|
2141 | div = ( u(k,j,i+1) - u(k,j,i) ) * ddx & |
---|
2142 | + ( v(k,j+1,i) - v(k,j,i) ) * ddy & |
---|
2143 | + ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) |
---|
2144 | |
---|
2145 | tend(k,j,i) = tend(k,j,i) - ( & |
---|
2146 | ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) - & |
---|
2147 | swap_diss_x_local(k,j) ) * ddx & |
---|
2148 | + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k) - & |
---|
2149 | swap_diss_y_local(k) ) * ddy & |
---|
2150 | + ( flux_t(k) + diss_t(k) - flux_d - diss_d & |
---|
2151 | ) * ddzw(k) & |
---|
2152 | ) + sk(k,j,i) * div |
---|
2153 | |
---|
2154 | swap_flux_y_local(k) = flux_n(k) |
---|
2155 | swap_diss_y_local(k) = diss_n(k) |
---|
2156 | swap_flux_x_local(k,j) = flux_r(k) |
---|
2157 | swap_diss_x_local(k,j) = diss_r(k) |
---|
2158 | flux_d = flux_t(k) |
---|
2159 | diss_d = diss_t(k) |
---|
2160 | |
---|
2161 | ENDDO |
---|
2162 | ! |
---|
2163 | !-- evaluation of statistics |
---|
2164 | SELECT CASE ( sk_char ) |
---|
2165 | |
---|
2166 | CASE ( 'pt' ) |
---|
2167 | DO k = nzb, nzt |
---|
2168 | sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) & |
---|
2169 | + ( flux_t(k) + diss_t(k) ) & |
---|
2170 | * weight_substep(intermediate_timestep_count) |
---|
2171 | ENDDO |
---|
2172 | CASE ( 'sa' ) |
---|
2173 | DO k = nzb, nzt |
---|
2174 | sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) & |
---|
2175 | + ( flux_t(k) + diss_t(k) ) & |
---|
2176 | * weight_substep(intermediate_timestep_count) |
---|
2177 | ENDDO |
---|
2178 | CASE ( 'q' ) |
---|
2179 | DO k = nzb, nzt |
---|
2180 | sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) & |
---|
2181 | + ( flux_t(k) + diss_t(k) ) & |
---|
2182 | * weight_substep(intermediate_timestep_count) |
---|
2183 | ENDDO |
---|
2184 | |
---|
2185 | END SELECT |
---|
2186 | |
---|
2187 | ENDDO |
---|
2188 | ENDDO |
---|
2189 | |
---|
2190 | END SUBROUTINE advec_s_ws |
---|
2191 | |
---|
2192 | |
---|
2193 | !------------------------------------------------------------------------------! |
---|
2194 | ! Advection of u - Call for all grid points |
---|
2195 | !------------------------------------------------------------------------------! |
---|
2196 | SUBROUTINE advec_u_ws |
---|
2197 | |
---|
2198 | USE arrays_3d |
---|
2199 | USE constants |
---|
2200 | USE control_parameters |
---|
2201 | USE grid_variables |
---|
2202 | USE indices |
---|
2203 | USE statistics |
---|
2204 | |
---|
2205 | IMPLICIT NONE |
---|
2206 | |
---|
2207 | INTEGER :: i, j, k, k_mm, k_pp, k_ppp, tn = 0 |
---|
2208 | REAL :: diss_d, div, flux_d, gu, gv, v_comp, w_comp |
---|
2209 | REAL, DIMENSION(nzb+1:nzt) :: swap_diss_y_local_u, swap_flux_y_local_u |
---|
2210 | REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local_u, & |
---|
2211 | swap_flux_x_local_u |
---|
2212 | REAL, DIMENSION(nzb:nzt) :: diss_n, diss_r, diss_t, flux_n, flux_r, & |
---|
2213 | flux_t, u_comp |
---|
2214 | |
---|
2215 | gu = 2.0 * u_gtrans |
---|
2216 | gv = 2.0 * v_gtrans |
---|
2217 | |
---|
2218 | ! |
---|
2219 | !-- Compute the fluxes for the whole left boundary of the processor domain. |
---|
2220 | i = nxlu |
---|
2221 | DO j = nys, nyn |
---|
2222 | DO k = nzb+1, nzb_max |
---|
2223 | |
---|
2224 | u_comp(k) = u(k,j,i) + u(k,j,i-1) - gu |
---|
2225 | swap_flux_x_local_u(k,j) = u_comp(k) * ( & |
---|
2226 | ( 37.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 & |
---|
2227 | + 7.0 * IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 & |
---|
2228 | + IBITS(wall_flags_0(k,j,i),9,1) * adv_mom_1 & |
---|
2229 | ) * & |
---|
2230 | ( u(k,j,i) + u(k,j,i-1) ) & |
---|
2231 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 & |
---|
2232 | + IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 & |
---|
2233 | ) * & |
---|
2234 | ( u(k,j,i+1) + u(k,j,i-2) ) & |
---|
2235 | + ( IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 & |
---|
2236 | ) * & |
---|
2237 | ( u(k,j,i+2) + u(k,j,i-3) ) & |
---|
2238 | ) |
---|
2239 | |
---|
2240 | swap_diss_x_local_u(k,j) = - ABS( u_comp(k) ) * ( & |
---|
2241 | ( 10.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 & |
---|
2242 | + 3.0 * IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 & |
---|
2243 | + IBITS(wall_flags_0(k,j,i),9,1) * adv_mom_1 & |
---|
2244 | ) * & |
---|
2245 | ( u(k,j,i) - u(k,j,i-1) ) & |
---|
2246 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 & |
---|
2247 | + IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 & |
---|
2248 | ) * & |
---|
2249 | ( u(k,j,i+1) - u(k,j,i-2) ) & |
---|
2250 | + ( IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 & |
---|
2251 | ) * & |
---|
2252 | ( u(k,j,i+2) - u(k,j,i-3) ) & |
---|
2253 | ) |
---|
2254 | |
---|
2255 | ENDDO |
---|
2256 | |
---|
2257 | DO k = nzb_max+1, nzt |
---|
2258 | |
---|
2259 | u_comp(k) = u(k,j,i) + u(k,j,i-1) - gu |
---|
2260 | swap_flux_x_local_u(k,j) = u_comp(k) * ( & |
---|
2261 | 37.0 * ( u(k,j,i) + u(k,j,i-1) ) & |
---|
2262 | - 8.0 * ( u(k,j,i+1) + u(k,j,i-2) ) & |
---|
2263 | + ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5 |
---|
2264 | swap_diss_x_local_u(k,j) = - ABS(u_comp(k)) * ( & |
---|
2265 | 10.0 * ( u(k,j,i) - u(k,j,i-1) ) & |
---|
2266 | - 5.0 * ( u(k,j,i+1) - u(k,j,i-2) ) & |
---|
2267 | + ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5 |
---|
2268 | |
---|
2269 | ENDDO |
---|
2270 | ENDDO |
---|
2271 | |
---|
2272 | DO i = nxlu, nxr |
---|
2273 | ! |
---|
2274 | !-- The following loop computes the fluxes for the south boundary points |
---|
2275 | j = nys |
---|
2276 | DO k = nzb+1, nzb_max |
---|
2277 | |
---|
2278 | v_comp = v(k,j,i) + v(k,j,i-1) - gv |
---|
2279 | swap_flux_y_local_u(k) = v_comp * ( & |
---|
2280 | ( 37.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 & |
---|
2281 | + 7.0 * IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 & |
---|
2282 | + IBITS(wall_flags_0(k,j,i),12,1) * adv_mom_1 & |
---|
2283 | ) * & |
---|
2284 | ( u(k,j,i) + u(k,j-1,i) ) & |
---|
2285 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 & |
---|
2286 | + IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 & |
---|
2287 | ) * & |
---|
2288 | ( u(k,j+1,i) + u(k,j-2,i) ) & |
---|
2289 | + ( IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 & |
---|
2290 | ) * & |
---|
2291 | ( u(k,j+2,i) + u(k,j-3,i) ) & |
---|
2292 | ) |
---|
2293 | |
---|
2294 | swap_diss_y_local_u(k) = - ABS ( v_comp ) * ( & |
---|
2295 | ( 10.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 & |
---|
2296 | + 3.0 * IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 & |
---|
2297 | + IBITS(wall_flags_0(k,j,i),12,1) * adv_mom_1 & |
---|
2298 | ) * & |
---|
2299 | ( u(k,j,i) - u(k,j-1,i) ) & |
---|
2300 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 & |
---|
2301 | + IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 & |
---|
2302 | ) * & |
---|
2303 | ( u(k,j+1,i) - u(k,j-2,i) ) & |
---|
2304 | + ( IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 & |
---|
2305 | ) * & |
---|
2306 | ( u(k,j+2,i) - u(k,j-3,i) ) & |
---|
2307 | ) |
---|
2308 | |
---|
2309 | ENDDO |
---|
2310 | |
---|
2311 | DO k = nzb_max+1, nzt |
---|
2312 | |
---|
2313 | v_comp = v(k,j,i) + v(k,j,i-1) - gv |
---|
2314 | swap_flux_y_local_u(k) = v_comp * ( & |
---|
2315 | 37.0 * ( u(k,j,i) + u(k,j-1,i) ) & |
---|
2316 | - 8.0 * ( u(k,j+1,i) + u(k,j-2,i) ) & |
---|
2317 | + ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5 |
---|
2318 | swap_diss_y_local_u(k) = - ABS(v_comp) * ( & |
---|
2319 | 10.0 * ( u(k,j,i) - u(k,j-1,i) ) & |
---|
2320 | - 5.0 * ( u(k,j+1,i) - u(k,j-2,i) ) & |
---|
2321 | + ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5 |
---|
2322 | |
---|
2323 | ENDDO |
---|
2324 | ! |
---|
2325 | !-- Computation of interior fluxes and tendency terms |
---|
2326 | DO j = nys, nyn |
---|
2327 | |
---|
2328 | flux_t(0) = 0.0 |
---|
2329 | diss_t(0) = 0.0 |
---|
2330 | flux_d = 0.0 |
---|
2331 | diss_d = 0.0 |
---|
2332 | |
---|
2333 | DO k = nzb+1, nzb_max |
---|
2334 | |
---|
2335 | u_comp(k) = u(k,j,i+1) + u(k,j,i) |
---|
2336 | flux_r(k) = ( u_comp(k) - gu ) * ( & |
---|
2337 | ( 37.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 & |
---|
2338 | + 7.0 * IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 & |
---|
2339 | + IBITS(wall_flags_0(k,j,i),9,1) * adv_mom_1 & |
---|
2340 | ) * & |
---|
2341 | ( u(k,j,i+1) + u(k,j,i) ) & |
---|
2342 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 & |
---|
2343 | + IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 & |
---|
2344 | ) * & |
---|
2345 | ( u(k,j,i+2) + u(k,j,i-1) ) & |
---|
2346 | + ( IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 & |
---|
2347 | ) * & |
---|
2348 | ( u(k,j,i+3) + u(k,j,i-2) ) & |
---|
2349 | ) |
---|
2350 | |
---|
2351 | diss_r(k) = - ABS( u_comp(k) - gu ) * ( & |
---|
2352 | ( 10.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 & |
---|
2353 | + 3.0 * IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 & |
---|
2354 | + IBITS(wall_flags_0(k,j,i),9,1) * adv_mom_1 & |
---|
2355 | ) * & |
---|
2356 | ( u(k,j,i+1) - u(k,j,i) ) & |
---|
2357 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 & |
---|
2358 | + IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 & |
---|
2359 | ) * & |
---|
2360 | ( u(k,j,i+2) - u(k,j,i-1) ) & |
---|
2361 | + ( IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 & |
---|
2362 | ) * & |
---|
2363 | ( u(k,j,i+3) - u(k,j,i-2) ) & |
---|
2364 | ) |
---|
2365 | |
---|
2366 | v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv |
---|
2367 | flux_n(k) = v_comp * ( & |
---|
2368 | ( 37.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 & |
---|
2369 | + 7.0 * IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 & |
---|
2370 | + IBITS(wall_flags_0(k,j,i),12,1) * adv_mom_1 & |
---|
2371 | ) * & |
---|
2372 | ( u(k,j+1,i) + u(k,j,i) ) & |
---|
2373 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 & |
---|
2374 | + IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 & |
---|
2375 | ) * & |
---|
2376 | ( u(k,j+2,i) + u(k,j-1,i) ) & |
---|
2377 | + ( IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 & |
---|
2378 | ) * & |
---|
2379 | ( u(k,j+3,i) + u(k,j-2,i) ) & |
---|
2380 | ) |
---|
2381 | |
---|
2382 | diss_n(k) = - ABS ( v_comp ) * ( & |
---|
2383 | ( 10.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 & |
---|
2384 | + 3.0 * IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 & |
---|
2385 | + IBITS(wall_flags_0(k,j,i),12,1) * adv_mom_1 & |
---|
2386 | ) * & |
---|
2387 | ( u(k,j+1,i) - u(k,j,i) ) & |
---|
2388 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 & |
---|
2389 | + IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 & |
---|
2390 | ) * & |
---|
2391 | ( u(k,j+2,i) - u(k,j-1,i) ) & |
---|
2392 | + ( IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 & |
---|
2393 | ) * & |
---|
2394 | ( u(k,j+3,i) - u(k,j-2,i) ) & |
---|
2395 | ) |
---|
2396 | ! |
---|
2397 | !-- k index has to be modified near bottom and top, else array |
---|
2398 | !-- subscripts will be exceeded. |
---|
2399 | k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),17,1) |
---|
2400 | k_pp = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),15,1) ) |
---|
2401 | k_mm = k - 2 * IBITS(wall_flags_0(k,j,i),17,1) |
---|
2402 | |
---|
2403 | w_comp = w(k,j,i) + w(k,j,i-1) |
---|
2404 | flux_t(k) = w_comp * ( & |
---|
2405 | ( 37.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 & |
---|
2406 | + 7.0 * IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3 & |
---|
2407 | + IBITS(wall_flags_0(k,j,i),15,1) * adv_mom_1 & |
---|
2408 | ) * & |
---|
2409 | ( u(k+1,j,i) + u(k,j,i) ) & |
---|
2410 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 & |
---|
2411 | + IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3 & |
---|
2412 | ) * & |
---|
2413 | ( u(k_pp,j,i) + u(k-1,j,i) ) & |
---|
2414 | + ( IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 & |
---|
2415 | ) * & |
---|
2416 | ( u(k_ppp,j,i) + u(k_mm,j,i) ) & |
---|
2417 | ) |
---|
2418 | |
---|
2419 | diss_t(k) = - ABS( w_comp ) * ( & |
---|
2420 | ( 10.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 & |
---|
2421 | + 3.0 * IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3 & |
---|
2422 | + IBITS(wall_flags_0(k,j,i),15,1) * adv_mom_1 & |
---|
2423 | ) * & |
---|
2424 | ( u(k+1,j,i) - u(k,j,i) ) & |
---|
2425 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 & |
---|
2426 | + IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3 & |
---|
2427 | ) * & |
---|
2428 | ( u(k_pp,j,i) - u(k-1,j,i) ) & |
---|
2429 | + ( IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 & |
---|
2430 | ) * & |
---|
2431 | ( u(k_ppp,j,i) - u(k_mm,j,i) ) & |
---|
2432 | ) |
---|
2433 | ! |
---|
2434 | !-- Calculate the divergence of the velocity field. A respective |
---|
2435 | !-- correction is needed to overcome numerical instabilities caused |
---|
2436 | !-- by a not sufficient reduction of divergences near topography. |
---|
2437 | div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx & |
---|
2438 | + ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy & |
---|
2439 | + ( w_comp - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) & |
---|
2440 | * ddzw(k) & |
---|
2441 | ) * 0.5 |
---|
2442 | |
---|
2443 | tend(k,j,i) = tend(k,j,i) - ( & |
---|
2444 | ( flux_r(k) + diss_r(k) & |
---|
2445 | - swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx & |
---|
2446 | + ( flux_n(k) + diss_n(k) & |
---|
2447 | - swap_flux_y_local_u(k) - swap_diss_y_local_u(k) ) * ddy & |
---|
2448 | + ( flux_t(k) + diss_t(k) & |
---|
2449 | - flux_d - diss_d & |
---|
2450 | ) * ddzw(k) & |
---|
2451 | ) + div * u(k,j,i) |
---|
2452 | |
---|
2453 | swap_flux_x_local_u(k,j) = flux_r(k) |
---|
2454 | swap_diss_x_local_u(k,j) = diss_r(k) |
---|
2455 | swap_flux_y_local_u(k) = flux_n(k) |
---|
2456 | swap_diss_y_local_u(k) = diss_n(k) |
---|
2457 | flux_d = flux_t(k) |
---|
2458 | diss_d = diss_t(k) |
---|
2459 | ! |
---|
2460 | !-- Statistical Evaluation of u'u'. The factor has to be applied |
---|
2461 | !-- for right evaluation when gallilei_trans = .T. . |
---|
2462 | sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) & |
---|
2463 | + ( flux_r(k) * & |
---|
2464 | ( u_comp(k) - 2.0 * hom(k,1,1,0) ) & |
---|
2465 | / ( u_comp(k) - gu + 1.0E-20 ) & |
---|
2466 | + diss_r(k) * & |
---|
2467 | ABS( u_comp(k) - 2.0 * hom(k,1,1,0) ) & |
---|
2468 | / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) ) & |
---|
2469 | * weight_substep(intermediate_timestep_count) |
---|
2470 | ! |
---|
2471 | !-- Statistical Evaluation of w'u'. |
---|
2472 | sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) & |
---|
2473 | + ( flux_t(k) + diss_t(k) ) & |
---|
2474 | * weight_substep(intermediate_timestep_count) |
---|
2475 | ENDDO |
---|
2476 | |
---|
2477 | DO k = nzb_max+1, nzt |
---|
2478 | |
---|
2479 | u_comp(k) = u(k,j,i+1) + u(k,j,i) |
---|
2480 | flux_r(k) = ( u_comp(k) - gu ) * ( & |
---|
2481 | 37.0 * ( u(k,j,i+1) + u(k,j,i) ) & |
---|
2482 | - 8.0 * ( u(k,j,i+2) + u(k,j,i-1) ) & |
---|
2483 | + ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5 |
---|
2484 | diss_r(k) = - ABS( u_comp(k) - gu ) * ( & |
---|
2485 | 10.0 * ( u(k,j,i+1) - u(k,j,i) ) & |
---|
2486 | - 5.0 * ( u(k,j,i+2) - u(k,j,i-1) ) & |
---|
2487 | + ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5 |
---|
2488 | |
---|
2489 | v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv |
---|
2490 | flux_n(k) = v_comp * ( & |
---|
2491 | 37.0 * ( u(k,j+1,i) + u(k,j,i) ) & |
---|
2492 | - 8.0 * ( u(k,j+2,i) + u(k,j-1,i) ) & |
---|
2493 | + ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5 |
---|
2494 | diss_n(k) = - ABS( v_comp ) * ( & |
---|
2495 | 10.0 * ( u(k,j+1,i) - u(k,j,i) ) & |
---|
2496 | - 5.0 * ( u(k,j+2,i) - u(k,j-1,i) ) & |
---|
2497 | + ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5 |
---|
2498 | ! |
---|
2499 | !-- k index has to be modified near bottom and top, else array |
---|
2500 | !-- subscripts will be exceeded. |
---|
2501 | k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),17,1) |
---|
2502 | k_pp = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),15,1) ) |
---|
2503 | k_mm = k - 2 * IBITS(wall_flags_0(k,j,i),17,1) |
---|
2504 | |
---|
2505 | w_comp = w(k,j,i) + w(k,j,i-1) |
---|
2506 | flux_t(k) = w_comp * ( & |
---|
2507 | ( 37.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 & |
---|
2508 | + 7.0 * IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3 & |
---|
2509 | + IBITS(wall_flags_0(k,j,i),15,1) * adv_mom_1 & |
---|
2510 | ) * & |
---|
2511 | ( u(k+1,j,i) + u(k,j,i) ) & |
---|
2512 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 & |
---|
2513 | + IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3 & |
---|
2514 | ) * & |
---|
2515 | ( u(k_pp,j,i) + u(k-1,j,i) ) & |
---|
2516 | + ( IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 & |
---|
2517 | ) * & |
---|
2518 | ( u(k_ppp,j,i) + u(k_mm,j,i) ) & |
---|
2519 | ) |
---|
2520 | |
---|
2521 | diss_t(k) = - ABS( w_comp ) * ( & |
---|
2522 | ( 10.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 & |
---|
2523 | + 3.0 * IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3 & |
---|
2524 | + IBITS(wall_flags_0(k,j,i),15,1) * adv_mom_1 & |
---|
2525 | ) * & |
---|
2526 | ( u(k+1,j,i) - u(k,j,i) ) & |
---|
2527 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 & |
---|
2528 | + IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3 & |
---|
2529 | ) * & |
---|
2530 | ( u(k_pp,j,i) - u(k-1,j,i) ) & |
---|
2531 | + ( IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 & |
---|
2532 | ) * & |
---|
2533 | ( u(k_ppp,j,i) - u(k_mm,j,i) ) & |
---|
2534 | ) |
---|
2535 | ! |
---|
2536 | !-- Calculate the divergence of the velocity field. A respective |
---|
2537 | !-- correction is needed to overcome numerical instabilities caused |
---|
2538 | !-- by a not sufficient reduction of divergences near topography. |
---|
2539 | div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx & |
---|
2540 | + ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy & |
---|
2541 | + ( w_comp - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) & |
---|
2542 | * ddzw(k) & |
---|
2543 | ) * 0.5 |
---|
2544 | |
---|
2545 | tend(k,j,i) = tend(k,j,i) - ( & |
---|
2546 | ( flux_r(k) + diss_r(k) & |
---|
2547 | - swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx & |
---|
2548 | + ( flux_n(k) + diss_n(k) & |
---|
2549 | - swap_flux_y_local_u(k) - swap_diss_y_local_u(k) ) * ddy & |
---|
2550 | + ( flux_t(k) + diss_t(k) & |
---|
2551 | - flux_d - diss_d & |
---|
2552 | ) * ddzw(k) & |
---|
2553 | ) + div * u(k,j,i) |
---|
2554 | |
---|
2555 | swap_flux_x_local_u(k,j) = flux_r(k) |
---|
2556 | swap_diss_x_local_u(k,j) = diss_r(k) |
---|
2557 | swap_flux_y_local_u(k) = flux_n(k) |
---|
2558 | swap_diss_y_local_u(k) = diss_n(k) |
---|
2559 | flux_d = flux_t(k) |
---|
2560 | diss_d = diss_t(k) |
---|
2561 | ! |
---|
2562 | !-- Statistical Evaluation of u'u'. The factor has to be applied |
---|
2563 | !-- for right evaluation when gallilei_trans = .T. . |
---|
2564 | sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) & |
---|
2565 | + ( flux_r(k) * & |
---|
2566 | ( u_comp(k) - 2.0 * hom(k,1,1,0) ) & |
---|
2567 | / ( u_comp(k) - gu + 1.0E-20 ) & |
---|
2568 | + diss_r(k) * & |
---|
2569 | ABS( u_comp(k) - 2.0 * hom(k,1,1,0) ) & |
---|
2570 | / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) ) & |
---|
2571 | * weight_substep(intermediate_timestep_count) |
---|
2572 | ! |
---|
2573 | !-- Statistical Evaluation of w'u'. |
---|
2574 | sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) & |
---|
2575 | + ( flux_t(k) + diss_t(k) ) & |
---|
2576 | * weight_substep(intermediate_timestep_count) |
---|
2577 | ENDDO |
---|
2578 | ENDDO |
---|
2579 | ENDDO |
---|
2580 | sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn) |
---|
2581 | |
---|
2582 | |
---|
2583 | END SUBROUTINE advec_u_ws |
---|
2584 | |
---|
2585 | |
---|
2586 | !------------------------------------------------------------------------------! |
---|
2587 | ! Advection of v - Call for all grid points |
---|
2588 | !------------------------------------------------------------------------------! |
---|
2589 | SUBROUTINE advec_v_ws |
---|
2590 | |
---|
2591 | USE arrays_3d |
---|
2592 | USE constants |
---|
2593 | USE control_parameters |
---|
2594 | USE grid_variables |
---|
2595 | USE indices |
---|
2596 | USE statistics |
---|
2597 | |
---|
2598 | IMPLICIT NONE |
---|
2599 | |
---|
2600 | |
---|
2601 | INTEGER :: i, j, k, k_mm, k_pp, k_ppp, tn = 0 |
---|
2602 | REAL :: diss_d, div, flux_d, gu, gv, u_comp, w_comp |
---|
2603 | REAL, DIMENSION(nzb+1:nzt) :: swap_diss_y_local_v, swap_flux_y_local_v |
---|
2604 | REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local_v, & |
---|
2605 | swap_flux_x_local_v |
---|
2606 | REAL, DIMENSION(nzb:nzt) :: diss_n, diss_r, diss_t, flux_n, flux_r, & |
---|
2607 | flux_t, v_comp |
---|
2608 | |
---|
2609 | gu = 2.0 * u_gtrans |
---|
2610 | gv = 2.0 * v_gtrans |
---|
2611 | ! |
---|
2612 | !-- First compute the whole left boundary of the processor domain |
---|
2613 | i = nxl |
---|
2614 | DO j = nysv, nyn |
---|
2615 | DO k = nzb+1, nzb_max |
---|
2616 | |
---|
2617 | u_comp = u(k,j-1,i) + u(k,j,i) - gu |
---|
2618 | swap_flux_x_local_v(k,j) = u_comp * ( & |
---|
2619 | ( 37.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 & |
---|
2620 | + 7.0 * IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3 & |
---|
2621 | + IBITS(wall_flags_0(k,j,i),18,1) * adv_mom_1 & |
---|
2622 | ) * & |
---|
2623 | ( v(k,j,i) + v(k,j,i-1) ) & |
---|
2624 | - ( 8.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 & |
---|
2625 | + IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3 & |
---|
2626 | ) * & |
---|
2627 | ( v(k,j,i+1) + v(k,j,i-2) ) & |
---|
2628 | + ( IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 & |
---|
2629 | ) * & |
---|
2630 | ( v(k,j,i+2) + v(k,j,i-3) ) & |
---|
2631 | ) |
---|
2632 | |
---|
2633 | swap_diss_x_local_v(k,j) = - ABS( u_comp ) * ( & |
---|
2634 | ( 10.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 & |
---|
2635 | + 3.0 * IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3 & |
---|
2636 | + IBITS(wall_flags_0(k,j,i),18,1) * adv_mom_1 & |
---|
2637 | ) * & |
---|
2638 | ( v(k,j,i) - v(k,j,i-1) ) & |
---|
2639 | - ( 5.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 & |
---|
2640 | + IBITS(wall_flags_0(k,j,i),19,1) * |
---|