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

Last change on this file since 274 was 274, checked in by heinze, 12 years ago

Indentation of the message calls corrected

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