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

Last change on this file since 673 was 673, checked in by suehring, 11 years ago

Right computation of the pressure using Runge-Kutta weighting coefficients. Consideration of the pressure gradient during the time integration removed. Removed bugfix concerning velocity variances.

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