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

Last change on this file since 624 was 624, checked in by heinze, 11 years ago

Calculation of q*2 added

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