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

Last change on this file since 110 was 110, checked in by raasch, 14 years ago

New:
---
Allows runs for a coupled atmosphere-ocean LES,
coupling frequency is controlled by new d3par-parameter dt_coupling,
the coupling mode (atmosphere_to_ocean or ocean_to_atmosphere) for the
respective processes is read from environment variable coupling_mode,
which is set by the mpiexec-command,
communication between the two models is done using the intercommunicator
comm_inter,
local files opened by the ocean model get the additional suffic "_O".
Assume saturation at k=nzb_s_inner(j,i) for atmosphere coupled to ocean.

A momentum flux can be set as top boundary condition using the new
inipar parameter top_momentumflux_u|v.

Non-cyclic boundary conditions can be used along all horizontal directions.

Quantities w*p* and w"e can be output as vertical profiles.

Initial profiles are reset to constant profiles in case that initializing_actions /= 'set_constant_profiles'. (init_rankine)

Optionally calculate km and kh from initial TKE e_init.

Changed:


Remaining variables iran changed to iran_part (advec_particles, init_particles).

In case that the presure solver is not called for every Runge-Kutta substep
(call_psolver_at_all_substeps = .F.), it is called after the first substep
instead of the last. In that case, random perturbations are also added to the
velocity field after the first substep.

Initialization of km,kh = 0.00001 for ocean = .T. (for ocean = .F. it remains 0.01).

Allow data_output_pr= q, wq, w"q", w*q* for humidity = .T. (instead of cloud_physics = .T.).

Errors:


Bugs from code parts for non-cyclic boundary conditions are removed: loops for
u and v are starting from index nxlu, nysv, respectively. The radiation boundary
condition is used for every Runge-Kutta substep. Velocity phase speeds for
the radiation boundary conditions are calculated for the first Runge-Kutta
substep only and reused for the further substeps. New arrays c_u, c_v, and c_w
are defined for this purpose. Several index errors are removed from the
radiation boundary condition code parts. Upper bounds for calculating
u_0 and v_0 (in production_e) are nxr+1 and nyn+1 because otherwise these
values are not available in case of non-cyclic boundary conditions.

+dots_num_palm in module user, +module netcdf_control in user_init (both in user_interface)

Bugfix: wrong sign removed from the buoyancy production term in the case use_reference = .T. (production_e)

Bugfix: Error message concerning output of particle concentration (pc) modified (check_parameters).

Bugfix: Rayleigh damping for ocean fixed.

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