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

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

New:
---

Changed:


PALM can be generally installed on any kind of Linux-, IBM-AIX-, or NEC-SX-system by adding appropriate settings to the configuration file.

Scripts are also running under the public domain ksh.

All system relevant compile and link options as well as the host identifier (local_host) are specified in the configuration file.

Filetransfer by ftp removed (options -f removed from mrun and mbuild).

Call of (system-)FLUSH routine moved to new routine local_flush.

return_addres and return_username are read from ENVPAR-NAMELIST-file instead of using local_getenv.

Preprocessor strings for different linux clusters changed to "lc", some preprocessor directives renamed (new: intel_openmp_bug), preprocessor directives for old systems removed

advec_particles, check_open, cpu_log, cpu_statistics, data_output_dvrp, flow_statistics, header, init_dvrp, init_particles, init_1d_model, init_dvrp, init_pegrid, local_getenv, local_system, local_tremain, local_tremain_ini, modules, palm, parin, run_control

new:
local_flush

mbuild, mrun

Errors:


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