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

Last change on this file since 102 was 102, checked in by raasch, 14 years ago

preliminary version for coupled runs

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