source: palm/tags/release-3.3/SOURCE/flow_statistics.f90 @ 99

Last change on this file since 99 was 98, checked in by raasch, 17 years ago

updating comments and rc-file

  • Property svn:keywords set to Id
File size: 41.8 KB
Line 
1 SUBROUTINE flow_statistics
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: flow_statistics.f90 98 2007-06-21 09:36:33Z 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
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             ENDDO
401
402             DO  k = nzb_diff_s_outer(j,i)-1, nzt_diff
403!
404!--             Heat flux w"pt"
405                sums_l(k,16,tn) = sums_l(k,16,tn)                              &
406                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
407                                               * ( pt(k+1,j,i) - pt(k,j,i) )   &
408                                               * ddzu(k+1) * rmask(j,i,sr)
409
410
411!
412!--             Salinity flux w"sa"
413                IF ( ocean )  THEN
414                   sums_l(k,65,tn) = sums_l(k,65,tn)                           &
415                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
416                                               * ( sa(k+1,j,i) - sa(k,j,i) )   &
417                                               * ddzu(k+1) * rmask(j,i,sr)
418                ENDIF
419
420!
421!--             Buoyancy flux, water flux (humidity flux) w"q"
422                IF ( humidity ) THEN
423                   sums_l(k,45,tn) = sums_l(k,45,tn)                           &
424                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
425                                               * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
426                                               * ddzu(k+1) * rmask(j,i,sr)
427                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
428                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
429                                               * ( q(k+1,j,i) - q(k,j,i) )     &
430                                               * ddzu(k+1) * rmask(j,i,sr)
431                   IF ( cloud_physics ) THEN
432                      sums_l(k,51,tn) = sums_l(k,51,tn)                        &
433                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
434                                               * ( ( q(k+1,j,i) - ql(k+1,j,i) )&
435                                                - ( q(k,j,i) - ql(k,j,i) ) )   &
436                                               * ddzu(k+1) * rmask(j,i,sr) 
437                   ENDIF
438                ENDIF
439
440!
441!--             Passive scalar flux
442                IF ( passive_scalar )  THEN
443                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
444                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
445                                               * ( q(k+1,j,i) - q(k,j,i) )     &
446                                               * ddzu(k+1) * rmask(j,i,sr)
447                ENDIF
448
449             ENDDO
450
451!
452!--          Subgridscale fluxes in the Prandtl layer
453             IF ( use_surface_fluxes )  THEN
454                sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
455                                    usws(j,i) * rmask(j,i,sr)     ! w"u"
456                sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
457                                    vsws(j,i) * rmask(j,i,sr)     ! w"v"
458                sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
459                                    shf(j,i)  * rmask(j,i,sr)     ! w"pt"
460                sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
461                                    0.0 * rmask(j,i,sr)           ! u"pt"
462                sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
463                                    0.0 * rmask(j,i,sr)           ! v"pt"
464                IF ( ocean )  THEN
465                   sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
466                                       saswsb(j,i) * rmask(j,i,sr)  ! w"sa"
467                ENDIF
468                IF ( humidity )  THEN
469                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
470                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
471                   IF ( cloud_physics )  THEN
472                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (           &
473                                          ( 1.0 + 0.61 * q(nzb,j,i) ) *   &
474                                          shf(j,i) + 0.61 * pt(nzb,j,i) * &
475                                                     qsws(j,i)            &
476                                                              )
477!
478!--                   Formula does not work if ql(nzb) /= 0.0
479                      sums_l(nzb,51,tn) = sums_l(nzb,51,tn) + &   ! w"q" (w"qv")
480                                          qsws(j,i) * rmask(j,i,sr)
481                   ENDIF
482                ENDIF
483                IF ( passive_scalar )  THEN
484                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
485                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
486                ENDIF
487             ENDIF
488
489!
490!--          Subgridscale fluxes at the top surface
491             IF ( use_top_fluxes )  THEN
492                sums_l(nzt,16,tn) = sums_l(nzt,16,tn) + &
493                                    tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
494                sums_l(nzt,58,tn) = sums_l(nzt,58,tn) + &
495                                    0.0 * rmask(j,i,sr)           ! u"pt"
496                sums_l(nzt,61,tn) = sums_l(nzt,61,tn) + &
497                                    0.0 * rmask(j,i,sr)           ! v"pt"
498                IF ( ocean )  THEN
499                   sums_l(nzt,65,tn) = sums_l(nzt,65,tn) + &
500                                       saswst(j,i) * rmask(j,i,sr)  ! w"sa"
501                ENDIF
502                IF ( humidity )  THEN
503                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
504                                       qswst(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
505                   IF ( cloud_physics )  THEN
506                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (           &
507                                          ( 1.0 + 0.61 * q(nzt,j,i) ) *   &
508                                          tswst(j,i) + 0.61 * pt(nzt,j,i) * &
509                                                     qsws(j,i)            &
510                                                              )
511!
512!--                   Formula does not work if ql(nzb) /= 0.0
513                      sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + &   ! w"q" (w"qv")
514                                          qswst(j,i) * rmask(j,i,sr)
515                   ENDIF
516                ENDIF
517                IF ( passive_scalar )  THEN
518                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
519                                       qswst(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
520                ENDIF
521             ENDIF
522
523!
524!--          Resolved fluxes (can be computed for all horizontal points)
525!--          NOTE: for simplicity, nzb_s_outer is used below, although strictly
526!--          ----  speaking the following k-loop would have to be split up and
527!--                rearranged according to the staggered grid.
528             DO  k = nzb_s_outer(j,i), nzt
529                ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
530                              u(k+1,j,i) - hom(k+1,1,1,sr) )
531                vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
532                              v(k+1,j,i) - hom(k+1,1,2,sr) )
533                pts = 0.5 * ( pt(k,j,i)   - hom(k,1,4,sr) + &
534                              pt(k+1,j,i) - hom(k+1,1,4,sr) )
535!
536!--             Momentum flux w*u*
537                sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5 *                     &
538                                                    ( w(k,j,i-1) + w(k,j,i) ) &
539                                                    * ust * rmask(j,i,sr)
540!
541!--             Momentum flux w*v*
542                sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5 *                     &
543                                                    ( w(k,j-1,i) + w(k,j,i) ) &
544                                                    * vst * rmask(j,i,sr)
545!
546!--             Heat flux w*pt*
547!--             The following formula (comment line, not executed) does not
548!--             work if applied to subregions
549!                sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5 *                     &
550!                                                    ( pt(k,j,i)+pt(k+1,j,i) ) &
551!                                                    * w(k,j,i) * rmask(j,i,sr)
552                sums_l(k,17,tn) = sums_l(k,17,tn) + pts * w(k,j,i) * &
553                                                    rmask(j,i,sr)
554!
555!--             Higher moments
556                sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 * &
557                                                    rmask(j,i,sr)
558                sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) * &
559                                                    rmask(j,i,sr)
560
561!
562!--             Salinity flux and density (density does not belong to here,
563!--             but so far there is no other suitable place to calculate)
564                IF ( ocean )  THEN
565                   pts = 0.5 * ( sa(k,j,i)   - hom(k,1,23,sr) + &
566                                 sa(k+1,j,i) - hom(k+1,1,23,sr) )
567                   sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) * &
568                                                       rmask(j,i,sr)
569                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) * &
570                                                       rmask(j,i,sr)
571                ENDIF
572
573!
574!--             Buoyancy flux, water flux, humidity flux and liquid water
575!--             content
576                IF ( humidity )  THEN
577                   pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) + &
578                                 vpt(k+1,j,i) - hom(k+1,1,44,sr) )
579                   sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
580                                                       rmask(j,i,sr)
581                   pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
582                                 q(k+1,j,i) - hom(k+1,1,41,sr) )
583                   sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
584                                                       rmask(j,i,sr)
585                   IF ( cloud_physics  .OR.  cloud_droplets )  THEN
586                      pts = 0.5 *                                            &
587                             ( ( q(k,j,i)   - ql(k,j,i)   ) - hom(k,1,42,sr) &
588                             + ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) )
589                      sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) * &
590                                                          rmask(j,i,sr)
591                      sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * &
592                                                          rmask(j,i,sr)
593                   ENDIF
594                ENDIF
595
596!
597!--             Passive scalar flux
598                IF ( passive_scalar )  THEN
599                   pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
600                                 q(k+1,j,i) - hom(k+1,1,41,sr) )
601                   sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
602                                                       rmask(j,i,sr)
603                ENDIF
604
605!
606!--             Energy flux w*e*
607                sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5 *           &
608                                              ( ust**2 + vst**2 + w(k,j,i)**2 )&
609                                              * rmask(j,i,sr)
610         
611             ENDDO
612          ENDDO
613       ENDDO
614
615!
616!--    Density at top follows Neumann condition
617       IF ( ocean )  sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
618
619!
620!--    Divergence of vertical flux of resolved scale energy and pressure
621!--    fluctuations. First calculate the products, then the divergence.
622!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
623       IF ( hom(nzb+1,2,55,0) /= 0.0 )  THEN
624
625          sums_ll = 0.0  ! local array
626
627          !$OMP DO
628          DO  i = nxl, nxr
629             DO  j = nys, nyn
630                DO  k = nzb_s_outer(j,i)+1, nzt
631
632                   sums_ll(k,1) = sums_ll(k,1) + 0.5 * w(k,j,i) * (        &
633                  ( 0.25 * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1)   &
634                           - 2.0 * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) )     &
635                           ) )**2                                          &
636                + ( 0.25 * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i)   &
637                           - 2.0 * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) )     &
638                           ) )**2                                          &
639                   + w(k,j,i)**2                                  )
640
641                   sums_ll(k,2) = sums_ll(k,2) + 0.5 * w(k,j,i) &
642                                               * ( p(k,j,i) + p(k+1,j,i) )
643
644                ENDDO
645             ENDDO
646          ENDDO
647          sums_ll(0,1)     = 0.0    ! because w is zero at the bottom
648          sums_ll(nzt+1,1) = 0.0
649          sums_ll(0,2)     = 0.0
650          sums_ll(nzt+1,2) = 0.0
651
652          DO  k = nzb_s_outer(j,i)+1, nzt
653             sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
654             sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
655          ENDDO
656          sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
657          sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
658
659       ENDIF
660
661!
662!--    Divergence of vertical flux of SGS TKE
663       IF ( hom(nzb+1,2,57,0) /= 0.0 )  THEN
664
665          !$OMP DO
666          DO  i = nxl, nxr
667             DO  j = nys, nyn
668                DO  k = nzb_s_outer(j,i)+1, nzt
669
670                   sums_l(k,57,tn) = sums_l(k,57,tn) + (                       &
671                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
672                 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
673                                                 ) * ddzw(k)
674
675                ENDDO
676             ENDDO
677          ENDDO
678          sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
679
680       ENDIF
681
682!
683!--    Horizontal heat fluxes (subgrid, resolved, total).
684!--    Do it only, if profiles shall be plotted.
685       IF ( hom(nzb+1,2,58,0) /= 0.0 ) THEN
686
687          !$OMP DO
688          DO  i = nxl, nxr
689             DO  j = nys, nyn
690                DO  k = nzb_s_outer(j,i)+1, nzt
691!
692!--                Subgrid horizontal heat fluxes u"pt", v"pt"
693                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5 *                   &
694                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
695                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
696                                                 * ddx * rmask(j,i,sr)
697                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5 *                   &
698                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
699                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
700                                                 * ddy * rmask(j,i,sr)
701!
702!--                Resolved horizontal heat fluxes u*pt*, v*pt*
703                   sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
704                                                  ( u(k,j,i) - hom(k,1,1,sr) ) &
705                                       * 0.5 * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
706                                                 pt(k,j,i)   - hom(k,1,4,sr) )
707                   pts = 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
708                                 pt(k,j,i)   - hom(k,1,4,sr) )
709                   sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
710                                                  ( v(k,j,i) - hom(k,1,2,sr) ) &
711                                       * 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
712                                                 pt(k,j,i)   - hom(k,1,4,sr) )
713                ENDDO
714             ENDDO
715          ENDDO
716!
717!--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
718          sums_l(nzb,58,tn) = 0.0
719          sums_l(nzb,59,tn) = 0.0
720          sums_l(nzb,60,tn) = 0.0
721          sums_l(nzb,61,tn) = 0.0
722          sums_l(nzb,62,tn) = 0.0
723          sums_l(nzb,63,tn) = 0.0
724
725       ENDIF
726
727!
728!--    Calculate the user-defined profiles
729       CALL user_statistics( 'profiles', sr, tn )
730       !$OMP END PARALLEL
731
732!
733!--    Summation of thread sums
734       IF ( threads_per_task > 1 )  THEN
735          DO  i = 1, threads_per_task-1
736             sums_l(:,3,0)          = sums_l(:,3,0) + sums_l(:,3,i)
737             sums_l(:,4:40,0)       = sums_l(:,4:40,0) + sums_l(:,4:40,i)
738             sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
739                                      sums_l(:,45:pr_palm,i)
740             IF ( max_pr_user > 0 )  THEN
741                sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
742                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
743                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
744             ENDIF
745          ENDDO
746       ENDIF
747
748#if defined( __parallel )
749!
750!--    Compute total sum from local sums
751       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
752                           MPI_SUM, comm2d, ierr )
753#else
754       sums = sums_l(:,:,0)
755#endif
756
757!
758!--    Final values are obtained by division by the total number of grid points
759!--    used for summation. After that store profiles.
760!--    Profiles:
761       DO  k = nzb, nzt+1
762          sums(k,:pr_palm-2)      = sums(k,:pr_palm-2) / ngp_2dh_outer(k,sr)
763       ENDDO
764!--    Upstream-parts
765       sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
766!--    u* and so on
767!--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
768!--    size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
769!--    above the topography, they are being divided by ngp_2dh(sr)
770       sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)    / &
771                                    ngp_2dh(sr)
772!--    eges, e*
773       sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
774                                    ngp_3d_inner(sr)
775!--    Old and new divergence
776       sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
777                                    ngp_3d_inner(sr)
778
779!--    User-defined profiles
780       IF ( max_pr_user > 0 )  THEN
781          DO  k = nzb, nzt+1
782             sums(k,pr_palm+1:pr_palm+max_pr_user) = &
783                                    sums(k,pr_palm+1:pr_palm+max_pr_user) / &
784                                    ngp_2dh_outer(k,sr)
785          ENDDO
786       ENDIF
787
788!
789!--    Collect horizontal average in hom.
790!--    Compute deduced averages (e.g. total heat flux)
791       hom(:,1,3,sr)  = sums(:,3)      ! w
792       hom(:,1,8,sr)  = sums(:,8)      ! e     profiles 5-7 are initial profiles
793       hom(:,1,9,sr)  = sums(:,9)      ! km
794       hom(:,1,10,sr) = sums(:,10)     ! kh
795       hom(:,1,11,sr) = sums(:,11)     ! l
796       hom(:,1,12,sr) = sums(:,12)     ! w"u"
797       hom(:,1,13,sr) = sums(:,13)     ! w*u*
798       hom(:,1,14,sr) = sums(:,14)     ! w"v"
799       hom(:,1,15,sr) = sums(:,15)     ! w*v*
800       hom(:,1,16,sr) = sums(:,16)     ! w"pt"
801       hom(:,1,17,sr) = sums(:,17)     ! w*pt*
802       hom(:,1,18,sr) = sums(:,16) + sums(:,17)    ! wpt
803       hom(:,1,19,sr) = sums(:,12) + sums(:,13)    ! wu
804       hom(:,1,20,sr) = sums(:,14) + sums(:,15)    ! wv
805       hom(:,1,21,sr) = sums(:,21)     ! w*pt*BC
806       hom(:,1,22,sr) = sums(:,16) + sums(:,21)    ! wptBC
807                                       ! profile 24 is initial profile (sa)
808                                       ! profiles 25-29 left empty for initial
809                                       ! profiles
810       hom(:,1,30,sr) = sums(:,30)     ! u*2
811       hom(:,1,31,sr) = sums(:,31)     ! v*2
812       hom(:,1,32,sr) = sums(:,32)     ! w*2
813       hom(:,1,33,sr) = sums(:,33)     ! pt*2
814       hom(:,1,34,sr) = sums(:,34)     ! e*
815       hom(:,1,35,sr) = sums(:,35)     ! w*2pt*
816       hom(:,1,36,sr) = sums(:,36)     ! w*pt*2
817       hom(:,1,37,sr) = sums(:,37)     ! w*e*
818       hom(:,1,38,sr) = sums(:,38)     ! w*3
819       hom(:,1,39,sr) = sums(:,38) / ( sums(:,32) + 1E-20 )**1.5    ! Sw
820       hom(:,1,40,sr) = sums(:,40)     ! p
821       hom(:,1,45,sr) = sums(:,45)     ! w"q"
822       hom(:,1,46,sr) = sums(:,46)     ! w*vpt*       
823       hom(:,1,47,sr) = sums(:,45) + sums(:,46)    ! wvpt
824       hom(:,1,48,sr) = sums(:,48)     ! w"q" (w"qv")
825       hom(:,1,49,sr) = sums(:,49)     ! w*q* (w*qv*)
826       hom(:,1,50,sr) = sums(:,48) + sums(:,49)    ! wq (wqv)
827       hom(:,1,51,sr) = sums(:,51)     ! w"qv"
828       hom(:,1,52,sr) = sums(:,52)     ! w*qv*       
829       hom(:,1,53,sr) = sums(:,52) + sums(:,51)    ! wq (wqv)
830       hom(:,1,54,sr) = sums(:,54)     ! ql
831       hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
832       hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
833       hom(:,1,57,sr) = sums(:,57)     ! w"e/dz
834       hom(:,1,58,sr) = sums(:,58)     ! u"pt"
835       hom(:,1,59,sr) = sums(:,59)     ! u*pt*
836       hom(:,1,60,sr) = sums(:,58) + sums(:,59)    ! upt_t
837       hom(:,1,61,sr) = sums(:,61)     ! v"pt"
838       hom(:,1,62,sr) = sums(:,62)     ! v*pt*
839       hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
840       hom(:,1,64,sr) = sums(:,64)     ! rho
841       hom(:,1,65,sr) = sums(:,65)     ! w"sa"
842       hom(:,1,66,sr) = sums(:,66)     ! w*sa*
843       hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
844
845       hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
846                                       ! upstream-parts u_x, u_y, u_z, v_x,
847                                       ! v_y, usw. (in last but one profile)
848       hom(:,1,pr_palm,sr) =   sums(:,pr_palm) 
849                                       ! u*, w'u', w'v', t* (in last profile)
850
851       IF ( max_pr_user > 0 )  THEN    ! user-defined profiles
852          hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
853                               sums(:,pr_palm+1:pr_palm+max_pr_user)
854       ENDIF
855
856!
857!--    Determine the boundary layer height using two different schemes.
858!--    First scheme: Starting from the Earth's (Ocean's) surface, look for the
859!--    first relative minimum (maximum) of the total heat flux.
860!--    The corresponding height is assumed as the boundary layer height, if it
861!--    is less than 1.5 times the height where the heat flux becomes negative
862!--    (positive) for the first time.
863!--    NOTE: This criterion is still capable of improving!
864       z_i(1) = 0.0
865       first = .TRUE.
866       IF ( ocean )  THEN
867          DO  k = nzt, nzb+1, -1
868             IF ( first .AND. hom(k,1,18,sr) < 0.0 )  THEN
869                first = .FALSE.
870                height = zw(k)
871             ENDIF
872             IF ( hom(k,1,18,sr) < 0.0  .AND. &
873                  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
874                IF ( zw(k) < 1.5 * height )  THEN
875                   z_i(1) = zw(k)
876                ELSE
877                   z_i(1) = height
878                ENDIF
879                EXIT
880             ENDIF
881          ENDDO
882       ELSE
883          DO  k = nzb, nzt-1
884             IF ( first .AND. hom(k,1,18,sr) < 0.0 )  THEN
885                first = .FALSE.
886                height = zw(k)
887             ENDIF
888             IF ( hom(k,1,18,sr) < 0.0  .AND. &
889                  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
890                IF ( zw(k) < 1.5 * height )  THEN
891                   z_i(1) = zw(k)
892                ELSE
893                   z_i(1) = height
894                ENDIF
895                EXIT
896             ENDIF
897          ENDDO
898       ENDIF
899
900!
901!--    Second scheme: Starting from the top/bottom model boundary, look for
902!--    the first characteristic kink in the temperature profile, where the
903!--    originally stable stratification notably weakens.
904       z_i(2) = 0.0
905       IF ( ocean )  THEN
906          DO  k = nzb+1, nzt-1
907             IF ( ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) > &
908                  2.0 * ( hom(k+1,1,4,sr) - hom(k,1,4,sr) ) )  THEN
909                z_i(2) = zu(k)
910                EXIT
911             ENDIF
912          ENDDO
913       ELSE
914          DO  k = nzt-1, nzb+1, -1
915             IF ( ( hom(k+1,1,4,sr) - hom(k,1,4,sr) ) > &
916                  2.0 * ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) )  THEN
917                z_i(2) = zu(k)
918                EXIT
919             ENDIF
920          ENDDO
921       ENDIF
922
923       hom(nzb+6,1,pr_palm,sr) = z_i(1)
924       hom(nzb+7,1,pr_palm,sr) = z_i(2)
925
926!
927!--    Computation of both the characteristic vertical velocity and
928!--    the characteristic convective boundary layer temperature.
929!--    The horizontal average at nzb+1 is input for the average temperature.
930       IF ( hom(nzb,1,18,sr) > 0.0  .AND.  z_i(1) /= 0.0 )  THEN
931          hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) * &
932                                       hom(nzb,1,18,sr) *      &
933                                       ABS( z_i(1) ) )**0.333333333
934!--       so far this only works if Prandtl layer is used
935          hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
936       ELSE
937          hom(nzb+8,1,pr_palm,sr)  = 0.0
938          hom(nzb+11,1,pr_palm,sr) = 0.0
939       ENDIF
940
941!
942!--    Collect the time series quantities
943       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)     ! E
944       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)     ! E*
945       ts_value(3,sr) = dt_3d
946       ts_value(4,sr) = hom(nzb,1,pr_palm,sr)       ! u*
947       ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)     ! th*
948       ts_value(6,sr) = u_max
949       ts_value(7,sr) = v_max
950       ts_value(8,sr) = w_max
951       ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)    ! new divergence
952       ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)    ! old Divergence
953       ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)    ! z_i(1)
954       ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)    ! z_i(2)
955       ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)    ! w*
956       ts_value(14,sr) = hom(nzb,1,16,sr)           ! w'pt'   at k=0
957       ts_value(15,sr) = hom(nzb+1,1,16,sr)         ! w'pt'   at k=1
958       ts_value(16,sr) = hom(nzb+1,1,18,sr)         ! wpt     at k=1
959       ts_value(17,sr) = hom(nzb,1,4,sr)            ! pt(0)
960       ts_value(18,sr) = hom(nzb+1,1,4,sr)          ! pt(zp)
961       ts_value(19,sr) = hom(nzb+9,1,pr_palm-1,sr)  ! splptx
962       ts_value(20,sr) = hom(nzb+10,1,pr_palm-1,sr) ! splpty
963       ts_value(21,sr) = hom(nzb+11,1,pr_palm-1,sr) ! splptz
964       IF ( ts_value(5,sr) /= 0.0 )  THEN
965          ts_value(22,sr) = ts_value(4,sr)**2 / &
966                            ( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
967       ELSE
968          ts_value(22,sr) = 10000.0
969       ENDIF
970
971!
972!--    Calculate additional statistics provided by the user interface
973       CALL user_statistics( 'time_series', sr, 0 )
974
975    ENDDO    ! loop of the subregions
976
977!
978!-- If required, sum up horizontal averages for subsequent time averaging
979    IF ( do_sum )  THEN
980       IF ( average_count_pr == 0 )  hom_sum = 0.0
981       hom_sum = hom_sum + hom(:,1,:,:)
982       average_count_pr = average_count_pr + 1
983       do_sum = .FALSE.
984    ENDIF
985
986!
987!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
988!-- This flag is reset after each time step in time_integration.
989    flow_statistics_called = .TRUE.
990
991    CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
992
993
994 END SUBROUTINE flow_statistics
995
996
997
Note: See TracBrowser for help on using the repository browser.