source: palm/trunk/SOURCE/flow_statistics.f90 @ 801

Last change on this file since 801 was 801, checked in by suehring, 12 years ago

Bugfix concerning OpenMP parallelization. Calculation of turbulent fluxes in advec_ws is now thread-safe.

  • Property svn:keywords set to Id
File size: 52.2 KB
Line 
1 SUBROUTINE flow_statistics
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6! Calculation of turbulent fluxes in advec_ws is now thread-safe.
7!
8! Former revisions:
9! -----------------
10! $Id: flow_statistics.f90 801 2012-01-10 17:30:36Z suehring $
11!
12! 743 2011-08-18 16:10:16Z suehring
13! Calculation of turbulent fluxes with WS-scheme only for the whole model
14! domain, not for user-defined subregions.
15!
16! 709 2011-03-30 09:31:40Z raasch
17! formatting adjustments
18!
19! 699 2011-03-22 17:52:22Z suehring
20! Bugfix in calculation of vertical velocity skewness. The added absolute value
21! avoid negative values in the root. Negative values of w'w' can occur at the
22! top or bottom of the model domain due to degrading the order of advection
23! scheme. Furthermore the calculation will be the same for all advection
24! schemes.
25!
26! 696 2011-03-18 07:03:49Z raasch
27! Bugfix: Summation of Wicker-Skamarock scheme fluxes and variances for all
28! threads
29!
30! 678 2011-02-02 14:31:56Z raasch
31! Bugfix in calculation of the divergence of vertical flux of resolved scale
32! energy, pressure fluctuations, and flux of pressure fluctuation itself
33!
34! 673 2011-01-18 16:19:48Z suehring
35! Top bc for the horizontal velocity variances added for ocean runs.
36! Setting the corresponding bottom bc moved to advec_ws.
37!
38! 667 2010-12-23 12:06:00Z suehring/gryschka
39! When advection is computed with ws-scheme, turbulent fluxes are already
40! computed in the respective advection routines and buffered in arrays
41! sums_xx_ws_l(). This is due to a consistent treatment of statistics with the
42! numerics and to avoid unphysical kinks near the surface.
43! So some if requests has to be done to dicern between fluxes from ws-scheme
44! other advection schemes.
45! Furthermore the computation of z_i is only done if the heat flux exceeds a
46! minimum value. This affects only simulations of a neutral boundary layer and
47! is due to reasons of computations in the advection scheme.
48!
49! 624 2010-12-10 11:46:30Z heinze
50! Calculation of q*2 added
51!
52! 622 2010-12-10 08:08:13Z raasch
53! optional barriers included in order to speed up collective operations
54!
55! 388 2009-09-23 09:40:33Z raasch
56! Vertical profiles of potential density and hydrostatic pressure are
57! calculated.
58! Added missing timeseries calculation of w"q"(0), moved timeseries q* to the
59! end.
60! Temperature gradient criterion for estimating the boundary layer height
61! replaced by the gradient criterion of Sullivan et al. (1998).
62! Output of messages replaced by message handling routine.
63!
64! 197 2008-09-16 15:29:03Z raasch
65! Spline timeseries splptx etc. removed, timeseries w'u', w'v', w'q' (k=0)
66! added,
67! bugfix: divide sums(k,8) (e) and sums(k,34) (e*) by ngp_2dh_s_inner(k,sr)
68! (like other scalars)
69!
70! 133 2007-11-20 10:10:53Z letzel
71! Vertical profiles now based on nzb_s_inner; they are divided by
72! ngp_2dh_s_inner (scalars, procucts of scalars) and ngp_2dh (staggered
73! velocity components and their products, procucts of scalars and velocity
74! components), respectively.
75!
76! 106 2007-08-16 14:30:26Z raasch
77! Prescribed momentum fluxes at the top surface are used,
78! profiles for w*p* and w"e are calculated
79!
80! 97 2007-06-21 08:23:15Z raasch
81! Statistics for ocean version (salinity, density) added,
82! calculation of z_i and Deardorff velocity scale adjusted to be used with
83! the ocean version
84!
85! 87 2007-05-22 15:46:47Z raasch
86! Two more arguments added to user_statistics, which is now also called for
87! user-defined profiles,
88! var_hom and var_sum renamed pr_palm
89!
90! 82 2007-04-16 15:40:52Z raasch
91! Cpp-directive lcmuk changed to intel_openmp_bug
92!
93! 75 2007-03-22 09:54:05Z raasch
94! Collection of time series quantities moved from routine flow_statistics to
95! here, routine user_statistics is called for each statistic region,
96! moisture renamed humidity
97!
98! 19 2007-02-23 04:53:48Z raasch
99! fluxes at top modified (tswst, qswst)
100!
101! RCS Log replace by Id keyword, revision history cleaned up
102!
103! Revision 1.41  2006/08/04 14:37:50  raasch
104! Error removed in non-parallel part (sums_l)
105!
106! Revision 1.1  1997/08/11 06:15:17  raasch
107! Initial revision
108!
109!
110! Description:
111! ------------
112! Compute average profiles and further average flow quantities for the different
113! user-defined (sub-)regions. The region indexed 0 is the total model domain.
114!
115! NOTE: For simplicity, nzb_s_inner and nzb_diff_s_inner are being used as a
116! ----  lower vertical index for k-loops for all variables, although strictly
117! speaking the k-loops would have to be split up according to the staggered
118! grid. However, this implies no error since staggered velocity components are
119! zero at the walls and inside buildings.
120!------------------------------------------------------------------------------!
121
122    USE arrays_3d
123    USE cloud_parameters
124    USE control_parameters
125    USE cpulog
126    USE grid_variables
127    USE indices
128    USE interfaces
129    USE pegrid
130    USE statistics
131
132    IMPLICIT NONE
133
134    INTEGER ::  i, j, k, omp_get_thread_num, sr, tn
135    LOGICAL ::  first
136    REAL    ::  dptdz_threshold, height, pts, sums_l_eper, sums_l_etot, ust, &
137                ust2, u2, vst, vst2, v2, w2, z_i(2)
138    REAL    ::  dptdz(nzb+1:nzt+1)
139    REAL    ::  sums_ll(nzb:nzt+1,2)
140
141    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
142
143!
144!-- To be on the safe side, check whether flow_statistics has already been
145!-- called once after the current time step
146    IF ( flow_statistics_called )  THEN
147
148       message_string = 'flow_statistics is called two times within one ' // &
149                        'timestep'
150       CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
151     
152    ENDIF
153
154!
155!-- Compute statistics for each (sub-)region
156    DO  sr = 0, statistic_regions
157
158!
159!--    Initialize (local) summation array
160       sums_l = 0.0
161
162!
163!--    Store sums that have been computed in other subroutines in summation
164!--    array
165       sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
166!--    WARNING: next line still has to be adjusted for OpenMP
167       sums_l(:,21,0) = sums_wsts_bc_l(:,sr)  ! heat flux from advec_s_bc
168       sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
169       sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
170
171!
172!--    Copy the turbulent quantities, evaluated in the advection routines to
173!--    the local array sums_l() for further computations
174       IF ( ws_scheme_mom .AND. sr == 0 )  THEN
175
176!       
177!--       According to the Neumann bc for the horizontal velocity components,
178!--       the corresponding fluxes has to satisfiy the same bc.
179          IF ( ocean )  THEN
180             sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
181             sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)   
182          ENDIF
183
184          DO  i = 0, threads_per_task-1
185!         
186!--          Swap the turbulent quantities evaluated in advec_ws.
187             sums_l(:,13,i) = sums_wsus_ws_l(:,i)       ! w*u*
188             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)       ! w*v*
189             sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2
190             sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
191             sums_l(:,32,i) = sums_ws2_ws_l(:,i)        ! w*2
192             sums_l(:,34,i) = sums_l(:,34,i) + 0.5 *                        & 
193                              ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +   &
194                                sums_ws2_ws_l(:,i) )    ! e*
195             DO  k = nzb, nzt
196                sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5 * ( &
197                                                      sums_us2_ws_l(k,i) +  &
198                                                      sums_vs2_ws_l(k,i) +  &
199                                                      sums_ws2_ws_l(k,i) )
200             ENDDO
201          ENDDO
202
203       ENDIF
204
205       IF ( ws_scheme_sca .AND. sr == 0 )  THEN
206
207          DO  i = 0, threads_per_task-1
208             sums_l(:,17,i) = sums_wspts_ws_l(:,i)      ! w*pt* from advec_s_ws
209             IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
210             IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) =              &
211                                                   sums_wsqs_ws_l(:,i) !w*q*
212          ENDDO
213
214       ENDIF
215!
216!--    Horizontally averaged profiles of horizontal velocities and temperature.
217!--    They must have been computed before, because they are already required
218!--    for other horizontal averages.
219       tn = 0
220
221       !$OMP PARALLEL PRIVATE( i, j, k, tn )
222#if defined( __intel_openmp_bug )
223       tn = omp_get_thread_num()
224#else
225!$     tn = omp_get_thread_num()
226#endif
227
228       !$OMP DO
229       DO  i = nxl, nxr
230          DO  j =  nys, nyn
231             DO  k = nzb_s_inner(j,i), nzt+1
232                sums_l(k,1,tn)  = sums_l(k,1,tn)  + u(k,j,i)  * rmask(j,i,sr)
233                sums_l(k,2,tn)  = sums_l(k,2,tn)  + v(k,j,i)  * rmask(j,i,sr)
234                sums_l(k,4,tn)  = sums_l(k,4,tn)  + pt(k,j,i) * rmask(j,i,sr)
235             ENDDO
236          ENDDO
237       ENDDO
238
239!
240!--    Horizontally averaged profile of salinity
241       IF ( ocean )  THEN
242          !$OMP DO
243          DO  i = nxl, nxr
244             DO  j =  nys, nyn
245                DO  k = nzb_s_inner(j,i), nzt+1
246                   sums_l(k,23,tn)  = sums_l(k,23,tn) + &
247                                      sa(k,j,i) * rmask(j,i,sr)
248                ENDDO
249             ENDDO
250          ENDDO
251       ENDIF
252
253!
254!--    Horizontally averaged profiles of virtual potential temperature,
255!--    total water content, specific humidity and liquid water potential
256!--    temperature
257       IF ( humidity )  THEN
258          !$OMP DO
259          DO  i = nxl, nxr
260             DO  j =  nys, nyn
261                DO  k = nzb_s_inner(j,i), nzt+1
262                   sums_l(k,44,tn)  = sums_l(k,44,tn) + &
263                                      vpt(k,j,i) * rmask(j,i,sr)
264                   sums_l(k,41,tn)  = sums_l(k,41,tn) + &
265                                      q(k,j,i) * rmask(j,i,sr)
266                ENDDO
267             ENDDO
268          ENDDO
269          IF ( cloud_physics )  THEN
270             !$OMP DO
271             DO  i = nxl, nxr
272                DO  j =  nys, nyn
273                   DO  k = nzb_s_inner(j,i), nzt+1
274                      sums_l(k,42,tn) = sums_l(k,42,tn) + &
275                                      ( q(k,j,i) - ql(k,j,i) ) * rmask(j,i,sr)
276                      sums_l(k,43,tn) = sums_l(k,43,tn) + ( &
277                                      pt(k,j,i) + l_d_cp*pt_d_t(k) * ql(k,j,i) &
278                                                          ) * rmask(j,i,sr)
279                   ENDDO
280                ENDDO
281             ENDDO
282          ENDIF
283       ENDIF
284
285!
286!--    Horizontally averaged profiles of passive scalar
287       IF ( passive_scalar )  THEN
288          !$OMP DO
289          DO  i = nxl, nxr
290             DO  j =  nys, nyn
291                DO  k = nzb_s_inner(j,i), nzt+1
292                   sums_l(k,41,tn)  = sums_l(k,41,tn) + q(k,j,i) * rmask(j,i,sr)
293                ENDDO
294             ENDDO
295          ENDDO
296       ENDIF
297       !$OMP END PARALLEL
298!
299!--    Summation of thread sums
300       IF ( threads_per_task > 1 )  THEN
301          DO  i = 1, threads_per_task-1
302             sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i)
303             sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
304             sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
305             IF ( ocean )  THEN
306                sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
307             ENDIF
308             IF ( humidity )  THEN
309                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
310                sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i)
311                IF ( cloud_physics )  THEN
312                   sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i)
313                   sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i)
314                ENDIF
315             ENDIF
316             IF ( passive_scalar )  THEN
317                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
318             ENDIF
319          ENDDO
320       ENDIF
321
322#if defined( __parallel )
323!
324!--    Compute total sum from local sums
325       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
326       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL, &
327                           MPI_SUM, comm2d, ierr )
328       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
329       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL, &
330                           MPI_SUM, comm2d, ierr )
331       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
332       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL, &
333                           MPI_SUM, comm2d, ierr )
334       IF ( ocean )  THEN
335          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
336          CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb, &
337                              MPI_REAL, MPI_SUM, comm2d, ierr )
338       ENDIF
339       IF ( humidity ) THEN
340          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
341          CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb, &
342                              MPI_REAL, MPI_SUM, comm2d, ierr )
343          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
344          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
345                              MPI_REAL, MPI_SUM, comm2d, ierr )
346          IF ( cloud_physics ) THEN
347             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
348             CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb, &
349                                 MPI_REAL, MPI_SUM, comm2d, ierr )
350             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
351             CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb, &
352                                 MPI_REAL, MPI_SUM, comm2d, ierr )
353          ENDIF
354       ENDIF
355
356       IF ( passive_scalar )  THEN
357          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
358          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
359                              MPI_REAL, MPI_SUM, comm2d, ierr )
360       ENDIF
361#else
362       sums(:,1) = sums_l(:,1,0)
363       sums(:,2) = sums_l(:,2,0)
364       sums(:,4) = sums_l(:,4,0)
365       IF ( ocean )  sums(:,23) = sums_l(:,23,0)
366       IF ( humidity ) THEN
367          sums(:,44) = sums_l(:,44,0)
368          sums(:,41) = sums_l(:,41,0)
369          IF ( cloud_physics ) THEN
370             sums(:,42) = sums_l(:,42,0)
371             sums(:,43) = sums_l(:,43,0)
372          ENDIF
373       ENDIF
374       IF ( passive_scalar )  sums(:,41) = sums_l(:,41,0)
375#endif
376
377!
378!--    Final values are obtained by division by the total number of grid points
379!--    used for summation. After that store profiles.
380       sums(:,1) = sums(:,1) / ngp_2dh(sr)
381       sums(:,2) = sums(:,2) / ngp_2dh(sr)
382       sums(:,4) = sums(:,4) / ngp_2dh_s_inner(:,sr)
383       hom(:,1,1,sr) = sums(:,1)             ! u
384       hom(:,1,2,sr) = sums(:,2)             ! v
385       hom(:,1,4,sr) = sums(:,4)             ! pt
386
387
388!
389!--    Salinity
390       IF ( ocean )  THEN
391          sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
392          hom(:,1,23,sr) = sums(:,23)             ! sa
393       ENDIF
394
395!
396!--    Humidity and cloud parameters
397       IF ( humidity ) THEN
398          sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
399          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
400          hom(:,1,44,sr) = sums(:,44)             ! vpt
401          hom(:,1,41,sr) = sums(:,41)             ! qv (q)
402          IF ( cloud_physics ) THEN
403             sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
404             sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
405             hom(:,1,42,sr) = sums(:,42)             ! qv
406             hom(:,1,43,sr) = sums(:,43)             ! pt
407          ENDIF
408       ENDIF
409
410!
411!--    Passive scalar
412       IF ( passive_scalar )  hom(:,1,41,sr) = sums(:,41) /  &
413            ngp_2dh_s_inner(:,sr)                    ! s (q)
414
415!
416!--    Horizontally averaged profiles of the remaining prognostic variables,
417!--    variances, the total and the perturbation energy (single values in last
418!--    column of sums_l) and some diagnostic quantities.
419!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
420!--    ----  speaking the following k-loop would have to be split up and
421!--          rearranged according to the staggered grid.
422!--          However, this implies no error since staggered velocity components
423!--          are zero at the walls and inside buildings.
424       tn = 0
425#if defined( __intel_openmp_bug )
426       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, &
427       !$OMP                    tn, ust, ust2, u2, vst, vst2, v2, w2 )
428       tn = omp_get_thread_num()
429#else
430       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2, w2 )
431!$     tn = omp_get_thread_num()
432#endif
433       !$OMP DO
434       DO  i = nxl, nxr
435          DO  j =  nys, nyn
436             sums_l_etot = 0.0
437             DO  k = nzb_s_inner(j,i), nzt+1
438!
439!--             Prognostic and diagnostic variables
440                sums_l(k,3,tn)  = sums_l(k,3,tn)  + w(k,j,i)  * rmask(j,i,sr)
441                sums_l(k,8,tn)  = sums_l(k,8,tn)  + e(k,j,i)  * rmask(j,i,sr)
442                sums_l(k,9,tn)  = sums_l(k,9,tn)  + km(k,j,i) * rmask(j,i,sr)
443                sums_l(k,10,tn) = sums_l(k,10,tn) + kh(k,j,i) * rmask(j,i,sr)
444                sums_l(k,40,tn) = sums_l(k,40,tn) + p(k,j,i)
445
446                sums_l(k,33,tn) = sums_l(k,33,tn) + &
447                                  ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr)
448
449                IF ( humidity )  THEN
450                   sums_l(k,70,tn) = sums_l(k,70,tn) + &
451                                  ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)
452                ENDIF
453               
454!
455!--             Higher moments
456!--             (Computation of the skewness of w further below)
457                sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i)**3 * rmask(j,i,sr)
458
459                sums_l_etot  = sums_l_etot + &
460                                        0.5 * ( u(k,j,i)**2 + v(k,j,i)**2 +    &
461                                        w(k,j,i)**2 ) * rmask(j,i,sr)
462             ENDDO
463!
464!--          Total and perturbation energy for the total domain (being
465!--          collected in the last column of sums_l). Summation of these
466!--          quantities is seperated from the previous loop in order to
467!--          allow vectorization of that loop.
468             sums_l(nzb+4,pr_palm,tn) = sums_l(nzb+4,pr_palm,tn) + sums_l_etot
469!
470!--          2D-arrays (being collected in the last column of sums_l)
471             sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +   &
472                                        us(j,i)   * rmask(j,i,sr)
473             sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) + &
474                                        usws(j,i) * rmask(j,i,sr)
475             sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) + &
476                                        vsws(j,i) * rmask(j,i,sr)
477             sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) + &
478                                        ts(j,i)   * rmask(j,i,sr)
479             IF ( humidity )  THEN
480                sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) + &
481                                            qs(j,i)   * rmask(j,i,sr)
482             ENDIF
483          ENDDO
484       ENDDO
485
486!
487!--    Computation of statistics when ws-scheme is not used. Else these
488!--    quantities are evaluated in the advection routines.
489       IF ( .NOT. ws_scheme_mom .OR. sr /= 0 )  THEN
490          !$OMP DO
491          DO  i = nxl, nxr
492             DO  j =  nys, nyn
493                sums_l_eper = 0.0
494                DO  k = nzb_s_inner(j,i), nzt+1
495                   u2   = u(k,j,i)**2
496                   v2   = v(k,j,i)**2
497                   w2   = w(k,j,i)**2
498                   ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
499                   vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
500
501                   sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)
502                   sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)
503                   sums_l(k,32,tn) = sums_l(k,32,tn) + w2   * rmask(j,i,sr)
504!
505!--             Perturbation energy
506
507                   sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5 *       &
508                                  ( ust2 + vst2 + w2 ) * rmask(j,i,sr)
509                   sums_l_eper  = sums_l_eper + &
510                                        0.5 * ( ust2+vst2+w2 ) * rmask(j,i,sr)
511
512                ENDDO
513                sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn)   &
514                     + sums_l_eper
515             ENDDO
516          ENDDO
517       ENDIF
518!
519!--    Horizontally averaged profiles of the vertical fluxes
520
521       !$OMP DO
522       DO  i = nxl, nxr
523          DO  j = nys, nyn
524!
525!--          Subgridscale fluxes (without Prandtl layer from k=nzb,
526!--          oterwise from k=nzb+1)
527!--          NOTE: for simplicity, nzb_diff_s_inner is used below, although
528!--          ----  strictly speaking the following k-loop would have to be
529!--                split up according to the staggered grid.
530!--                However, this implies no error since staggered velocity
531!--                components are zero at the walls and inside buildings.
532
533             DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
534!
535!--             Momentum flux w"u"
536                sums_l(k,12,tn) = sums_l(k,12,tn) - 0.25 * (                   &
537                               km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
538                                                           ) * (               &
539                                   ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
540                                 + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
541                                                               ) * rmask(j,i,sr)
542!
543!--             Momentum flux w"v"
544                sums_l(k,14,tn) = sums_l(k,14,tn) - 0.25 * (                   &
545                               km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
546                                                           ) * (               &
547                                   ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
548                                 + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
549                                                               ) * rmask(j,i,sr)
550!
551!--             Heat flux w"pt"
552                sums_l(k,16,tn) = sums_l(k,16,tn)                              &
553                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
554                                               * ( pt(k+1,j,i) - pt(k,j,i) )   &
555                                               * ddzu(k+1) * rmask(j,i,sr)
556
557
558!
559!--             Salinity flux w"sa"
560                IF ( ocean )  THEN
561                   sums_l(k,65,tn) = sums_l(k,65,tn)                           &
562                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
563                                               * ( sa(k+1,j,i) - sa(k,j,i) )   &
564                                               * ddzu(k+1) * rmask(j,i,sr)
565                ENDIF
566
567!
568!--             Buoyancy flux, water flux (humidity flux) w"q"
569                IF ( humidity ) THEN
570                   sums_l(k,45,tn) = sums_l(k,45,tn)                           &
571                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
572                                               * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
573                                               * ddzu(k+1) * rmask(j,i,sr)
574                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
575                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
576                                               * ( q(k+1,j,i) - q(k,j,i) )     &
577                                               * ddzu(k+1) * rmask(j,i,sr)
578                   IF ( cloud_physics ) THEN
579                      sums_l(k,51,tn) = sums_l(k,51,tn)                        &
580                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
581                                               * ( ( q(k+1,j,i) - ql(k+1,j,i) )&
582                                                - ( q(k,j,i) - ql(k,j,i) ) )   &
583                                               * ddzu(k+1) * rmask(j,i,sr) 
584                   ENDIF
585                ENDIF
586
587!
588!--             Passive scalar flux
589                IF ( passive_scalar )  THEN
590                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
591                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
592                                               * ( q(k+1,j,i) - q(k,j,i) )     &
593                                               * ddzu(k+1) * rmask(j,i,sr)
594                ENDIF
595
596             ENDDO
597
598!
599!--          Subgridscale fluxes in the Prandtl layer
600             IF ( use_surface_fluxes )  THEN
601                sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
602                                    usws(j,i) * rmask(j,i,sr)     ! w"u"
603                sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
604                                    vsws(j,i) * rmask(j,i,sr)     ! w"v"
605                sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
606                                    shf(j,i)  * rmask(j,i,sr)     ! w"pt"
607                sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
608                                    0.0 * rmask(j,i,sr)           ! u"pt"
609                sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
610                                    0.0 * rmask(j,i,sr)           ! v"pt"
611                IF ( ocean )  THEN
612                   sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
613                                       saswsb(j,i) * rmask(j,i,sr)  ! w"sa"
614                ENDIF
615                IF ( humidity )  THEN
616                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
617                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
618                   IF ( cloud_physics )  THEN
619                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (           &
620                                          ( 1.0 + 0.61 * q(nzb,j,i) ) *   &
621                                          shf(j,i) + 0.61 * pt(nzb,j,i) * &
622                                                     qsws(j,i)            &
623                                                              )
624!
625!--                   Formula does not work if ql(nzb) /= 0.0
626                      sums_l(nzb,51,tn) = sums_l(nzb,51,tn) + &   ! w"q" (w"qv")
627                                          qsws(j,i) * rmask(j,i,sr)
628                   ENDIF
629                ENDIF
630                IF ( passive_scalar )  THEN
631                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
632                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
633                ENDIF
634             ENDIF
635
636!
637!--          Subgridscale fluxes at the top surface
638             IF ( use_top_fluxes )  THEN
639                sums_l(nzt:nzt+1,12,tn) = sums_l(nzt:nzt+1,12,tn) + &
640                                    uswst(j,i) * rmask(j,i,sr)    ! w"u"
641                sums_l(nzt:nzt+1,14,tn) = sums_l(nzt:nzt+1,14,tn) + &
642                                    vswst(j,i) * rmask(j,i,sr)    ! w"v"
643                sums_l(nzt:nzt+1,16,tn) = sums_l(nzt:nzt+1,16,tn) + &
644                                    tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
645                sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
646                                    0.0 * rmask(j,i,sr)           ! u"pt"
647                sums_l(nzt:nzt+1,61,tn) = sums_l(nzt:nzt+1,61,tn) + &
648                                    0.0 * rmask(j,i,sr)           ! v"pt"
649
650                IF ( ocean )  THEN
651                   sums_l(nzt,65,tn) = sums_l(nzt,65,tn) + &
652                                       saswst(j,i) * rmask(j,i,sr)  ! w"sa"
653                ENDIF
654                IF ( humidity )  THEN
655                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
656                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
657                   IF ( cloud_physics )  THEN
658                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (           &
659                                          ( 1.0 + 0.61 * q(nzt,j,i) ) *   &
660                                          tswst(j,i) + 0.61 * pt(nzt,j,i) * &
661                                                     qsws(j,i)            &
662                                                              )
663!
664!--                   Formula does not work if ql(nzb) /= 0.0
665                      sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + &   ! w"q" (w"qv")
666                                          qswst(j,i) * rmask(j,i,sr)
667                   ENDIF
668                ENDIF
669                IF ( passive_scalar )  THEN
670                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
671                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
672                ENDIF
673             ENDIF
674
675!
676!--          Resolved fluxes (can be computed for all horizontal points)
677!--          NOTE: for simplicity, nzb_s_inner is used below, although strictly
678!--          ----  speaking the following k-loop would have to be split up and
679!--                rearranged according to the staggered grid.
680             DO  k = nzb_s_inner(j,i), nzt
681                ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
682                              u(k+1,j,i) - hom(k+1,1,1,sr) )
683                vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
684                              v(k+1,j,i) - hom(k+1,1,2,sr) )
685                pts = 0.5 * ( pt(k,j,i)   - hom(k,1,4,sr) + &
686                              pt(k+1,j,i) - hom(k+1,1,4,sr) )
687
688!--             Higher moments
689                sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 * &
690                                                    rmask(j,i,sr)
691                sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) * &
692                                                    rmask(j,i,sr)
693
694!
695!--             Salinity flux and density (density does not belong to here,
696!--             but so far there is no other suitable place to calculate)
697                IF ( ocean )  THEN
698                   IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
699                      pts = 0.5 * ( sa(k,j,i)   - hom(k,1,23,sr) + &
700                                 sa(k+1,j,i) - hom(k+1,1,23,sr) )
701                      sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) * &
702                                                       rmask(j,i,sr)
703                   ENDIF
704                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) * &
705                                                       rmask(j,i,sr)
706                   sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) * &
707                                                       rmask(j,i,sr)
708                ENDIF
709
710!
711!--             Buoyancy flux, water flux, humidity flux and liquid water
712!--             content
713                IF ( humidity )  THEN
714                   pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) + &
715                                 vpt(k+1,j,i) - hom(k+1,1,44,sr) )
716                   sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
717                                                       rmask(j,i,sr)
718
719                   IF ( cloud_physics  .OR.  cloud_droplets )  THEN
720                      pts = 0.5 *                                            &
721                             ( ( q(k,j,i)   - ql(k,j,i)   ) - hom(k,1,42,sr) &
722                             + ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) )
723                      sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) * &
724                                                          rmask(j,i,sr)
725                      sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * &
726                                                          rmask(j,i,sr)
727                   ENDIF
728                ENDIF
729
730!
731!--             Passive scalar flux
732                IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca   &
733                     .OR. sr /= 0 ) )  THEN
734                   pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
735                                 q(k+1,j,i) - hom(k+1,1,41,sr) )
736                   sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
737                                                       rmask(j,i,sr)
738                ENDIF
739
740!
741!--             Energy flux w*e*
742!--             has to be adjusted
743                sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5 *          &
744                                             ( ust**2 + vst**2 + w(k,j,i)**2 )&
745                                             * rmask(j,i,sr)
746             ENDDO
747          ENDDO
748       ENDDO
749!
750!--    For speed optimization fluxes which have been computed in part directly
751!--    inside the WS advection routines are treated seperatly
752!--    Momentum fluxes first:
753       IF ( .NOT. ws_scheme_mom .OR. sr /= 0  )  THEN
754         !$OMP DO
755         DO  i = nxl, nxr
756            DO  j = nys, nyn
757               DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
758                  ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
759                              u(k+1,j,i) - hom(k+1,1,1,sr) )
760                  vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
761                              v(k+1,j,i) - hom(k+1,1,2,sr) )
762!                             
763!--               Momentum flux w*u*
764                  sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5 *                   &
765                                                    ( w(k,j,i-1) + w(k,j,i) ) &
766                                                    * ust * rmask(j,i,sr)
767!
768!--               Momentum flux w*v*
769                  sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5 *                   &
770                                                    ( w(k,j-1,i) + w(k,j,i) ) &
771                                                    * vst * rmask(j,i,sr)
772               ENDDO
773            ENDDO
774         ENDDO
775
776       ENDIF
777       IF ( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
778         !$OMP DO
779         DO  i = nxl, nxr
780            DO  j = nys, nyn
781               DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
782!
783!--               Vertical heat flux
784                  sums_l(k,17,tn) = sums_l(k,17,tn) +  0.5 * &
785                           ( pt(k,j,i)   - hom(k,1,4,sr) + &
786                           pt(k+1,j,i) - hom(k+1,1,4,sr) ) &
787                           * w(k,j,i) * rmask(j,i,sr)
788                  IF ( humidity )  THEN
789                     pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
790                                q(k+1,j,i) - hom(k+1,1,41,sr) )
791                     sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
792                                                      rmask(j,i,sr)
793                  ENDIF
794               ENDDO
795            ENDDO
796         ENDDO
797
798       ENDIF
799
800
801!
802!--    Density at top follows Neumann condition
803       IF ( ocean )  THEN
804          sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
805          sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
806       ENDIF
807
808!
809!--    Divergence of vertical flux of resolved scale energy and pressure
810!--    fluctuations as well as flux of pressure fluctuation itself (68).
811!--    First calculate the products, then the divergence.
812!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
813       IF ( hom(nzb+1,2,55,0) /= 0.0  .OR.  hom(nzb+1,2,68,0) /= 0.0 )  THEN
814
815          sums_ll = 0.0  ! local array
816
817          !$OMP DO
818          DO  i = nxl, nxr
819             DO  j = nys, nyn
820                DO  k = nzb_s_inner(j,i)+1, nzt
821
822                   sums_ll(k,1) = sums_ll(k,1) + 0.5 * w(k,j,i) * (        &
823                  ( 0.25 * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1)   &
824                           - 0.5 * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) )     &
825                           ) )**2                                          &
826                + ( 0.25 * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i)   &
827                           - 0.5 * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) )     &
828                           ) )**2                                          &
829                   + w(k,j,i)**2                                  )
830
831                   sums_ll(k,2) = sums_ll(k,2) + 0.5 * w(k,j,i) &
832                                               * ( p(k,j,i) + p(k+1,j,i) )
833
834                ENDDO
835             ENDDO
836          ENDDO
837          sums_ll(0,1)     = 0.0    ! because w is zero at the bottom
838          sums_ll(nzt+1,1) = 0.0
839          sums_ll(0,2)     = 0.0
840          sums_ll(nzt+1,2) = 0.0
841
842          DO  k = nzb+1, nzt
843             sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
844             sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
845             sums_l(k,68,tn) = sums_ll(k,2)
846          ENDDO
847          sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
848          sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
849          sums_l(nzb,68,tn) = 0.0    ! because w* = 0 at nzb
850
851       ENDIF
852
853!
854!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
855       IF ( hom(nzb+1,2,57,0) /= 0.0  .OR.  hom(nzb+1,2,69,0) /= 0.0 )  THEN
856
857          !$OMP DO
858          DO  i = nxl, nxr
859             DO  j = nys, nyn
860                DO  k = nzb_s_inner(j,i)+1, nzt
861
862                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5 * (                 &
863                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
864                 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
865                                                             ) * ddzw(k)
866
867                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5 * (                 &
868                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
869                                                              )
870
871                ENDDO
872             ENDDO
873          ENDDO
874          sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
875          sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
876
877       ENDIF
878
879!
880!--    Horizontal heat fluxes (subgrid, resolved, total).
881!--    Do it only, if profiles shall be plotted.
882       IF ( hom(nzb+1,2,58,0) /= 0.0 ) THEN
883
884          !$OMP DO
885          DO  i = nxl, nxr
886             DO  j = nys, nyn
887                DO  k = nzb_s_inner(j,i)+1, nzt
888!
889!--                Subgrid horizontal heat fluxes u"pt", v"pt"
890                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5 *                   &
891                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
892                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
893                                                 * ddx * rmask(j,i,sr)
894                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5 *                   &
895                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
896                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
897                                                 * ddy * rmask(j,i,sr)
898!
899!--                Resolved horizontal heat fluxes u*pt*, v*pt*
900                   sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
901                                                  ( u(k,j,i) - hom(k,1,1,sr) ) &
902                                       * 0.5 * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
903                                                 pt(k,j,i)   - hom(k,1,4,sr) )
904                   pts = 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
905                                 pt(k,j,i)   - hom(k,1,4,sr) )
906                   sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
907                                                  ( v(k,j,i) - hom(k,1,2,sr) ) &
908                                       * 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
909                                                 pt(k,j,i)   - hom(k,1,4,sr) )
910                ENDDO
911             ENDDO
912          ENDDO
913!
914!--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
915          sums_l(nzb,58,tn) = 0.0
916          sums_l(nzb,59,tn) = 0.0
917          sums_l(nzb,60,tn) = 0.0
918          sums_l(nzb,61,tn) = 0.0
919          sums_l(nzb,62,tn) = 0.0
920          sums_l(nzb,63,tn) = 0.0
921
922       ENDIF
923
924!
925!--    Calculate the user-defined profiles
926       CALL user_statistics( 'profiles', sr, tn )
927       !$OMP END PARALLEL
928
929!
930!--    Summation of thread sums
931       IF ( threads_per_task > 1 )  THEN
932          DO  i = 1, threads_per_task-1
933             sums_l(:,3,0)          = sums_l(:,3,0) + sums_l(:,3,i)
934             sums_l(:,4:40,0)       = sums_l(:,4:40,0) + sums_l(:,4:40,i)
935             sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
936                                      sums_l(:,45:pr_palm,i)
937             IF ( max_pr_user > 0 )  THEN
938                sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
939                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
940                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
941             ENDIF
942          ENDDO
943       ENDIF
944
945#if defined( __parallel )
946
947!
948!--    Compute total sum from local sums
949       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
950       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
951                           MPI_SUM, comm2d, ierr )
952#else
953       sums = sums_l(:,:,0)
954#endif
955
956!
957!--    Final values are obtained by division by the total number of grid points
958!--    used for summation. After that store profiles.
959!--    Profiles:
960       DO  k = nzb, nzt+1
961          sums(k,3)               = sums(k,3)           / ngp_2dh(sr)
962          sums(k,8:11)            = sums(k,8:11)        / ngp_2dh_s_inner(k,sr)
963          sums(k,12:22)           = sums(k,12:22)       / ngp_2dh(sr)
964          sums(k,23:29)           = sums(k,23:29)       / ngp_2dh_s_inner(k,sr)
965          sums(k,30:32)           = sums(k,30:32)       / ngp_2dh(sr)
966          sums(k,33:34)           = sums(k,33:34)       / ngp_2dh_s_inner(k,sr)
967          sums(k,35:39)           = sums(k,35:39)       / ngp_2dh(sr)
968          sums(k,40)              = sums(k,40)          / ngp_2dh_s_inner(k,sr)
969          sums(k,45:53)           = sums(k,45:53)       / ngp_2dh(sr)
970          sums(k,54)              = sums(k,54)          / ngp_2dh_s_inner(k,sr)
971          sums(k,55:63)           = sums(k,55:63)       / ngp_2dh(sr)
972          sums(k,64)              = sums(k,64)          / ngp_2dh_s_inner(k,sr)
973          sums(k,65:69)           = sums(k,65:69)       / ngp_2dh(sr)
974          sums(k,70:pr_palm-2)    = sums(k,70:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
975       ENDDO
976
977!--    Upstream-parts
978       sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
979!--    u* and so on
980!--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
981!--    size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
982!--    above the topography, they are being divided by ngp_2dh(sr)
983       sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)    / &
984                                    ngp_2dh(sr)
985       sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
986                                    ngp_2dh(sr)
987!--    eges, e*
988       sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
989                                    ngp_3d(sr)
990!--    Old and new divergence
991       sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
992                                    ngp_3d_inner(sr)
993
994!--    User-defined profiles
995       IF ( max_pr_user > 0 )  THEN
996          DO  k = nzb, nzt+1
997             sums(k,pr_palm+1:pr_palm+max_pr_user) = &
998                                    sums(k,pr_palm+1:pr_palm+max_pr_user) / &
999                                    ngp_2dh_s_inner(k,sr)
1000          ENDDO
1001       ENDIF
1002!
1003!--    Collect horizontal average in hom.
1004!--    Compute deduced averages (e.g. total heat flux)
1005       hom(:,1,3,sr)  = sums(:,3)      ! w
1006       hom(:,1,8,sr)  = sums(:,8)      ! e     profiles 5-7 are initial profiles
1007       hom(:,1,9,sr)  = sums(:,9)      ! km
1008       hom(:,1,10,sr) = sums(:,10)     ! kh
1009       hom(:,1,11,sr) = sums(:,11)     ! l
1010       hom(:,1,12,sr) = sums(:,12)     ! w"u"
1011       hom(:,1,13,sr) = sums(:,13)     ! w*u*
1012       hom(:,1,14,sr) = sums(:,14)     ! w"v"
1013       hom(:,1,15,sr) = sums(:,15)     ! w*v*
1014       hom(:,1,16,sr) = sums(:,16)     ! w"pt"
1015       hom(:,1,17,sr) = sums(:,17)     ! w*pt*
1016       hom(:,1,18,sr) = sums(:,16) + sums(:,17)    ! wpt
1017       hom(:,1,19,sr) = sums(:,12) + sums(:,13)    ! wu
1018       hom(:,1,20,sr) = sums(:,14) + sums(:,15)    ! wv
1019       hom(:,1,21,sr) = sums(:,21)     ! w*pt*BC
1020       hom(:,1,22,sr) = sums(:,16) + sums(:,21)    ! wptBC
1021                                       ! profile 24 is initial profile (sa)
1022                                       ! profiles 25-29 left empty for initial
1023                                       ! profiles
1024       hom(:,1,30,sr) = sums(:,30)     ! u*2
1025       hom(:,1,31,sr) = sums(:,31)     ! v*2
1026       hom(:,1,32,sr) = sums(:,32)     ! w*2
1027       hom(:,1,33,sr) = sums(:,33)     ! pt*2
1028       hom(:,1,34,sr) = sums(:,34)     ! e*
1029       hom(:,1,35,sr) = sums(:,35)     ! w*2pt*
1030       hom(:,1,36,sr) = sums(:,36)     ! w*pt*2
1031       hom(:,1,37,sr) = sums(:,37)     ! w*e*
1032       hom(:,1,38,sr) = sums(:,38)     ! w*3
1033       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20 )**1.5   ! Sw
1034       hom(:,1,40,sr) = sums(:,40)     ! p
1035       hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
1036       hom(:,1,46,sr) = sums(:,46)     ! w*vpt*       
1037       hom(:,1,47,sr) = sums(:,45) + sums(:,46)    ! wvpt
1038       hom(:,1,48,sr) = sums(:,48)     ! w"q" (w"qv")
1039       hom(:,1,49,sr) = sums(:,49)     ! w*q* (w*qv*)
1040       hom(:,1,50,sr) = sums(:,48) + sums(:,49)    ! wq (wqv)
1041       hom(:,1,51,sr) = sums(:,51)     ! w"qv"
1042       hom(:,1,52,sr) = sums(:,52)     ! w*qv*       
1043       hom(:,1,53,sr) = sums(:,52) + sums(:,51)    ! wq (wqv)
1044       hom(:,1,54,sr) = sums(:,54)     ! ql
1045       hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
1046       hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
1047       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho )/dz
1048       hom(:,1,58,sr) = sums(:,58)     ! u"pt"
1049       hom(:,1,59,sr) = sums(:,59)     ! u*pt*
1050       hom(:,1,60,sr) = sums(:,58) + sums(:,59)    ! upt_t
1051       hom(:,1,61,sr) = sums(:,61)     ! v"pt"
1052       hom(:,1,62,sr) = sums(:,62)     ! v*pt*
1053       hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
1054       hom(:,1,64,sr) = sums(:,64)     ! rho
1055       hom(:,1,65,sr) = sums(:,65)     ! w"sa"
1056       hom(:,1,66,sr) = sums(:,66)     ! w*sa*
1057       hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
1058       hom(:,1,68,sr) = sums(:,68)     ! w*p*
1059       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho
1060       hom(:,1,70,sr) = sums(:,70)     ! q*2
1061       hom(:,1,71,sr) = sums(:,71)     ! prho
1062       hom(:,1,72,sr) = hyp * 1E-4     ! hyp in dbar
1063
1064       hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
1065                                       ! upstream-parts u_x, u_y, u_z, v_x,
1066                                       ! v_y, usw. (in last but one profile)
1067       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
1068                                       ! u*, w'u', w'v', t* (in last profile)
1069
1070       IF ( max_pr_user > 0 )  THEN    ! user-defined profiles
1071          hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
1072                               sums(:,pr_palm+1:pr_palm+max_pr_user)
1073       ENDIF
1074
1075!
1076!--    Determine the boundary layer height using two different schemes.
1077!--    First scheme: Starting from the Earth's (Ocean's) surface, look for the
1078!--    first relative minimum (maximum) of the total heat flux.
1079!--    The corresponding height is assumed as the boundary layer height, if it
1080!--    is less than 1.5 times the height where the heat flux becomes negative
1081!--    (positive) for the first time.
1082       z_i(1) = 0.0
1083       first = .TRUE.
1084
1085       IF ( ocean )  THEN
1086          DO  k = nzt, nzb+1, -1
1087             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
1088                .AND. abs(hom(k,1,18,sr)) > 1.0E-8)  THEN
1089                first = .FALSE.
1090                height = zw(k)
1091             ENDIF
1092             IF ( hom(k,1,18,sr) < 0.0  .AND. &
1093                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
1094                  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
1095                IF ( zw(k) < 1.5 * height )  THEN
1096                   z_i(1) = zw(k)
1097                ELSE
1098                   z_i(1) = height
1099                ENDIF
1100                EXIT
1101             ENDIF
1102          ENDDO
1103       ELSE
1104          DO  k = nzb, nzt-1
1105             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
1106                .AND. abs(hom(k,1,18,sr)) > 1.0E-8 )  THEN
1107                first = .FALSE.
1108                height = zw(k)
1109             ENDIF
1110             IF ( hom(k,1,18,sr) < 0.0  .AND. & 
1111                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
1112                  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
1113                IF ( zw(k) < 1.5 * height )  THEN
1114                   z_i(1) = zw(k)
1115                ELSE
1116                   z_i(1) = height
1117                ENDIF
1118                EXIT
1119             ENDIF
1120          ENDDO
1121       ENDIF
1122
1123!
1124!--    Second scheme: Gradient scheme from Sullivan et al. (1998), modified
1125!--    by Uhlenbrock(2006). The boundary layer height is the height with the
1126!--    maximal local temperature gradient: starting from the second (the last
1127!--    but one) vertical gridpoint, the local gradient must be at least
1128!--    0.2K/100m and greater than the next four gradients.
1129!--    WARNING: The threshold value of 0.2K/100m must be adjusted for the
1130!--             ocean case!
1131       z_i(2) = 0.0
1132       DO  k = nzb+1, nzt+1
1133          dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
1134       ENDDO
1135       dptdz_threshold = 0.2 / 100.0
1136
1137       IF ( ocean )  THEN
1138          DO  k = nzt+1, nzb+5, -1
1139             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
1140                  dptdz(k) > dptdz(k-1)  .AND.  dptdz(k) > dptdz(k-2)  .AND.  &
1141                  dptdz(k) > dptdz(k-3)  .AND.  dptdz(k) > dptdz(k-4) )  THEN
1142                z_i(2) = zw(k-1)
1143                EXIT
1144             ENDIF
1145          ENDDO
1146       ELSE
1147          DO  k = nzb+1, nzt-3
1148             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
1149                  dptdz(k) > dptdz(k+1)  .AND.  dptdz(k) > dptdz(k+2)  .AND.  &
1150                  dptdz(k) > dptdz(k+3)  .AND.  dptdz(k) > dptdz(k+4) )  THEN
1151                z_i(2) = zw(k-1)
1152                EXIT
1153             ENDIF
1154          ENDDO
1155       ENDIF
1156
1157       hom(nzb+6,1,pr_palm,sr) = z_i(1)
1158       hom(nzb+7,1,pr_palm,sr) = z_i(2)
1159
1160!
1161!--    Computation of both the characteristic vertical velocity and
1162!--    the characteristic convective boundary layer temperature.
1163!--    The horizontal average at nzb+1 is input for the average temperature.
1164       IF ( hom(nzb,1,18,sr) > 0.0 .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8 &
1165           .AND.  z_i(1) /= 0.0 )  THEN
1166          hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) * &
1167                                       hom(nzb,1,18,sr) *      &
1168                                       ABS( z_i(1) ) )**0.333333333
1169!--       so far this only works if Prandtl layer is used
1170          hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
1171       ELSE
1172          hom(nzb+8,1,pr_palm,sr)  = 0.0
1173          hom(nzb+11,1,pr_palm,sr) = 0.0
1174       ENDIF
1175
1176!
1177!--    Collect the time series quantities
1178       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)     ! E
1179       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)     ! E*
1180       ts_value(3,sr) = dt_3d
1181       ts_value(4,sr) = hom(nzb,1,pr_palm,sr)       ! u*
1182       ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)     ! th*
1183       ts_value(6,sr) = u_max
1184       ts_value(7,sr) = v_max
1185       ts_value(8,sr) = w_max
1186       ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)    ! new divergence
1187       ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)    ! old Divergence
1188       ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)    ! z_i(1)
1189       ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)    ! z_i(2)
1190       ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)    ! w*
1191       ts_value(14,sr) = hom(nzb,1,16,sr)           ! w'pt'   at k=0
1192       ts_value(15,sr) = hom(nzb+1,1,16,sr)         ! w'pt'   at k=1
1193       ts_value(16,sr) = hom(nzb+1,1,18,sr)         ! wpt     at k=1
1194       ts_value(17,sr) = hom(nzb,1,4,sr)            ! pt(0)
1195       ts_value(18,sr) = hom(nzb+1,1,4,sr)          ! pt(zp)
1196       ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr)    ! u'w'    at k=0
1197       ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)    ! v'w'    at k=0
1198       ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
1199
1200       IF ( ts_value(5,sr) /= 0.0 )  THEN
1201          ts_value(22,sr) = ts_value(4,sr)**2 / &
1202                            ( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
1203       ELSE
1204          ts_value(22,sr) = 10000.0
1205       ENDIF
1206
1207       ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
1208!
1209!--    Calculate additional statistics provided by the user interface
1210       CALL user_statistics( 'time_series', sr, 0 )
1211
1212    ENDDO    ! loop of the subregions
1213
1214!
1215!-- If required, sum up horizontal averages for subsequent time averaging
1216    IF ( do_sum )  THEN
1217       IF ( average_count_pr == 0 )  hom_sum = 0.0
1218       hom_sum = hom_sum + hom(:,1,:,:)
1219       average_count_pr = average_count_pr + 1
1220       do_sum = .FALSE.
1221    ENDIF
1222
1223!
1224!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
1225!-- This flag is reset after each time step in time_integration.
1226    flow_statistics_called = .TRUE.
1227
1228    CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
1229
1230
1231 END SUBROUTINE flow_statistics
Note: See TracBrowser for help on using the repository browser.