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

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

last commit documented

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