source: palm/tags/release-3.2b/SOURCE/flow_statistics.f90 @ 139

Last change on this file since 139 was 90, checked in by raasch, 17 years ago

New:
---
Calculation and output of user-defined profiles. New &userpar parameters data_output_pr_user and max_pr_user.

check_parameters, flow_statistics, modules, parin, read_var_list, user_interface, write_var_list

Changed:


Division through dt_3d replaced by multiplication of the inverse. For performance optimisation, this is done in the loop calculating the divergence instead of using a seperate loop. (pres.f90) var_hom and var_sum renamed pr_palm.

data_output_profiles, flow_statistics, init_3d_model, modules, parin, pres, read_var_list, run_control, time_integration

Errors:


Bugfix: work_fft*_vec removed from some PRIVATE-declarations (poisfft).

Bugfix: field_chr renamed field_char (user_interface).

Bugfix: output of use_upstream_for_tke (header).

header, poisfft, user_interface

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