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

Last change on this file since 197 was 197, checked in by raasch, 13 years ago

further adjustments for SGI and other small changes

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